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Using a diagrammatic approach to Eulerian perturbation theory, we analytically calculate the 
variance and skewness of the density and velocity divergence induced by gravitational evolution 
from Gaussian initial conditions, including corrections beyond leading order. Except for the power 
spectrum, previous calculations in cosmological perturbation theory have been confined to leading 
order (tree level) — we extend these to include loop corrections. For scale-free initial power spectra, 
P(k) ~ k n with —2 < n < 2, the one-loop variance a 2 = (S 2 ) = a 2 + 1.82aj and the skewness 
S 3 = {6 3 }/(T 4 = 34/7 + 9.8(7?, where ai is the rms fluctuation of the density field to linear order. 
(These results depend weakly on the spectral index n, due to the non-locality of the non- linear 
solutions to the equations of motion.) Thus, loop corrections for the (unsmoothed) density field 
begin to dominate over tree-level contributions (and perturbation theory presumably begins to 
break down) when a 2 ~ 1/2. For the divergence of the velocity field, loop dominance does not 
occur until a 2 PS 1. We also compute loop corrections to the variance, skewness, and kurtosis for 
several non-linear approximation schemes, where the calculation can be easily generalized to 1-point 
cumulants of higher order and arbitrary number of loops. We find that the Zel'dovich approximation 
gives the best approximation to the loop corrections of exact perturbation theory, followed by the 
Linear Potential approximation (LPA) and the Frozen Flow approximation (FFA), in qualitative 
agreement with the relative behavior of tree-level results. In LPA and FFA, loop corrections are 
infrared divergent for spectral indices n < — 1; this is related to the breaking of Galilean invariance 
in these schemes. 



I. INTRODUCTION 

There is growing evidence that the large-scale structure of the Universe grew via gravitational instability from 
small primordial fluctuations in the matter density. At early epochs, the growth of density perturbations can be 
described by linear perturbation theory, provided that the linear power spectrum P{k) falls off" less steeply than k 4 
for small k [1-3]. In the linear regime, perturbation Fourier modes evolve independently of one another, conserving 
the statistical properties of the primordial fluctuations. In particular, if the primordial fluctuations are Gaussian 
random fields (as expected from the simplest inflationary models), they remain Gaussian in linear theory. In this 
case, the statistical properties of the density and velocity fields are completely determined by the two-point correlation 
function or the power spectrum. When the fluctuations become non-linear, coupling between different Fourier modes 
becomes important, inducing non-trivial correlations that modify the statistical properties of the cosmological fields. 
For Gaussian initial conditions, this causes the appearance of higher-order reduced correlations (p-point cumulants), 
which constitute independent statistics that can be measured in observational data and numerical simulations, even 
when the departure from the linear regime is small. 

Non-linear cosmological perturbation theory (NLCPT) provides a framework for analytic calculations of these 
effects in this weakly non-linear regime. The assumption is that even when the density and velocity fields are highly 
non-linear on the smallest scales, one can accurately describe the evolution of the large-scale properties of these fields 
during the first stages of non-linear evolution using a systematic perturbative approach. The validity of this ansatz 
has not yet been proven theoretically, but its predictions have been checked against numerical simulations. In fact, 
comparison of one-point statistics with N-body simulations shows that leading order NLCPT provides an adequate 
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description even on scales where next-to-leading order and higher order perturbative contributions would be expected 
to become important. The surprising agreement of leading order calculations with numerical simulations raises 
interesting questions: does it imply that the next-to-leading order corrections are generically small, or are they merely 
'accidentally' small for the cases that have been numerically studied? Or does it instead imply that the whole series 
of contributions beyond leading order conspire to nearly cancel? The answer to this puzzle should tell us something 
interesting about gravitational dynamics and the transition to the non-linear regime. These questions clearly motivate 
the study of corrections beyond leading order (loop corrections) in NLCPT. In any case, the calculation of next-to- 
leading order corrections should help us to better understand the successes and limitations of the perturbative approach 
to gravitational instability. 

To study the dynamics of gravitational instability, one must choose between an Eulerian or Lagrangian description; 
each has its advantages and drawbacks. In the Eulerian approach, one considers the motion of matter relative to some 
fixed coordinate system, whereas in the Lagrangian formalism one follows the trajectories of individual fluid elements. 
Perturbation theory is formulated quite differently in the two approaches: in the Eulerian framework, the small 
expansion parameter is the linear rms density fluctuation; in the Lagrangian formulation, the expansion is done about 
'inertial' fluid element trajectories. As a result, the Lagrangian description is more successful in describing large density 
contrasts that can arise before multi-streaming (orbit crossing) occurs; for an account of Lagrangian perturbation 
theory, see [4] and references therein. On the other hand, the Eulerian approach is conceptually straightforward and 
has the advantage that it works directly with the density and velocity fields whose statistics we want to quantify. 
In this paper, we use the Eulerian description, which can be solved to arbitrary order in perturbation theory for an 
Einstein-de Sitter Universe [5]. 

Leading order calculations in non-linear theory began with Peebles [6], who used second order perturbation theory to 
obtain the skewness for the unsmoothed density field (assuming Gaussian initial conditions), S3 = {& 3 ) I {& 2 ) 2 = 34/7. 
Fry [7] used second and third order perturbation theory to calculate the leading order three and four-point functions, 
respectively. He also showed how one can associate tree diagrams with each perturbative contribution to leading order. 
The predictions for the three-point function were subsequently found to agree with numerical simulations [8]. Goroff 
et al. [5] developed a systematic approach to NLCPT, obtaining recursion relations for the density and velocity fields 
to arbitrary order in an Einstein-de Sitter universe. They showed how the perturbation series could be represented in 
terms of Feynman diagrams, with leading order contributions corresponding to tree diagrams, next-to-leading order 
to one-loop diagrams, and so on. They also introduced the normalized one-point cumulants S p (p > 3), which are 
measures of the skewness (p = 3), kurtosis (p = 4), and higher central moments of the density field. They obtained 
the leading-order S p for p = 3, 4, 5 for the Gaussian-smoothed density field by Monte Carlo integration for SI = 1 cold 
dark matter (CDM), but they did not do any explicit loop calculations. These calculations were extended to non-flat 
Universes (£1 ^ 1) using Lagrangian perturbation theory [9,10], with the result that the skewness of the unsmoothed 
density field depends very weakly on SI. Recently, equivalent results were obtained in Eulerian perturbation theory [11]. 

To compare the results of NLCPT with numerical simulations and with observations, the one-point cumulants must 
be calculated for the density field smoothed over a finite volume. Analytic leading-order results for S3 were obtained 
in [12] for scale-free initial power spectra, P(k) ~ k n , for Gaussian (n > —1) and top-hat smoothing (—3 < n < 0); 
with smoothing, S3 decreases with increasing n. Comparison of these results with N-body simulations showed a very 
good agreement at a 2 « 0.5, while Monte Carlo integration of the smoothed skewness for standard CDM agreed with 
N-body results up to scales where a 2 k, 1. These results were later generalized by Bernardeau [13], who presented 
analytic calculations of S3, S^, T3 and T4 (T p being the normalized cumulants for the divergence of the velocity field) 
for top-hat smoothing using Eulerian perturbation theory for arbitrary SI, A, and P(k). The dependence on A was 
found to be extremely small, with only T3 and T4 being sensitive to SI (see also [14]). For Gaussian smoothing, S4 
was obtained by Monte Carlo integration in [15], and semianalytically in [16], where similar results are presented for 
T4 and also analytic computations of S3 and T3. 

An important generalization of all these results is the calculation of the whole series (i.e., for all p) of one-point 
cumulants at tree-level. This was done for unsmoothed fields by Bernardeau [17,18] and generalized to include top- 
hat smoothing for S p and T p [19], who showed that the tree-level cumulant generating function obeys the equations 
of motion of the spherical collapse model. Therefore, all the information concerning one-point cumulants at tree- 
level is encoded in the spherical collapse dynamics. In particular, the tree-level smoothed one-point cumulants are 
independent of the variance cr(R) and depend on the smoothing scale R only through derivatives of the spectral index 
n with respect to R; thus, for scale-free linear power spectra, the leading-order S P (R) are constants. From the series 
of leading-order cumulants, the one-point probability function P(S)dS was reconstructed and shown to agree with 
CDM N-body simulations even when a « 1.5 [20]. A detailed comparison between leading-order NLCPT predictions 
and N-body simulations for the S p parameters (up to p = 10) was carried out in [21,22], showing very good agreement 
up to scales where a « 1. 

One-loop corrections in NLCPT have been studied in the literature for the 2-point cumulant in real space (the 
two-point correlation function) and Fourier space (the power spectrum) [23-30]. These calculations show how power 
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is transferred between large and small scales; the agreement with N-body simulations is good for a 2 < 1 (when 
long wavelength modes dominate the non-linear contribution) and gradually breaks down when a increases above 
unity [29,30]. Higher order loop corrections to the power spectrum were considered in [31], including the full contri- 
butions up to 2 loops and the most important terms for high k in 3- and 4-loop contributions. Very recently, one-loop 
corrections to the variance were studied numerically in [32] for gaussian smoothing. 

Comparison of NLCPT predictions with galaxy surveys is more intricate. Due to the possibility of a non-trivial 
bias mechanism during the process of structure formation, the statistical properties of the galaxy distribution are not 
necessarily equivalent to those of the underlying matter distribution predicted by theory. Moreover, when considering 
higher-order correlations, the relation between the galaxy and matter density correlations is complicated by the 
fact that the bias between galaxies and mass is likely to be non-linear [33]. However, assuming the validity of 
NLCPT (as demonstrated by comparison with numerical simulations), one can use the observed galaxy higher- 
order correlations [34-37] to put constraints on the non-linear bias [38,39] and test models of large-scale structure 
formation [40,41]. Another possibility is to test perturbation theory by including some form of bias in numerical 
simulations [42]. 

Calculation of loop corrections to one-point cumulants constitutes a natural development of the subject and should 
provide further insight into the weakly non-linear regime. Numerical simulations have shown that the S p parameters 
depart from the tree-level perturbation theory predictions in the non-linear regime [43-45]. This is to be expected, 
since the fluid equations of motion which are the starting point for the perturbative expansion become invalid when 
multistreaming and vorticity develop. However, it is important to know whether NLCPT can provide an understanding 
of this departure in the first stages of non-linear evolution or whether this behavior is instead due to non-perturbative 
and/or physical effects not included in the analytic treatments so far. 

In this paper, we develop the diagrammatic approach to Eulerian perturbation theory, providing the calculational 
rules and technical machinery to perform loop calculations. As applications, we calculate analytically the one-loop 
corrections to the variance and skewness for the unsmoothed density and divergence of the velocity fields, assuming 
Gaussian initial conditions. Since we deal with unsmoothed fields, these results cannot yet be checked against 
numerical simulations. In subsequent papers, we will present results for the one-loop bispectrum (the three-point 
function in Fourier space) and the variance and skewness for the smoothed density field. 

When fluctuations become strongly non-linear, perturbation theory breaks down and one must rely on N-body 
simulations to follow the dynamical evolution. Alternatively, one can formulate analytic approximation schemes 
that provide some insight into the physics of the non-linear regime. This was first suggested by Zel'dovich [46], 
who initiated the Lagrangian approach to gravitational instability by extrapolating the linear solution in Lagrangian 
space into the non-linear regime. The Zel'dovich approximation gives an analytic understanding of the formation 
of walls and filaments in gravitational clustering. Recently, alternative approximations have been suggested, e.g., to 
extrapolate the linear velocity field into the non-linear regime [47] (the frozen flow approximation), and to extrapolate 
the linear gravitational potential [48,49] (the linear potential approximation). This set of non-linear approximations 
have in common that they are all consistent in linear theory with the exact dynamics, their differences appearing at 
second order in a perturbative expansion. 

We also calculate loop corrections to one-point cumulants in these non-linear approximation schemes; these calcu- 
lations are much easier to perform than for the exact dynamics and can be carried out to higher order and number of 
loops. Although these approximations do not provide a very accurate estimation of the exact one-point cumulants (as 
was first noticed for the Zel'dovich approximation [50]), they nevertheless give order of magnitude estimates which are 
useful when corresponding results for the exact dynamics are unknown. A comparison of the tree-level predictions in 
these approximations for one-point cumulants has been given in [51-53]. A comparative study of non-linear approxi- 
mations in the strongly non-linear regime was carried out in [54]. One-loop corrections to S3 and S4 for the Zel'dovich 
approximation in the unsmoothed case where estimated in [20] by taking moments of the (regularized) one-point 
probability distribution. Here we provide a perturbative calculation and extend the results to the divergence of the 
velocity field and to two-loop corrections for S3 and T3. We also calculate up to 3-loop corrections to the variance of 
the density and divergence of the velocity fields. 

The paper is organized as follows. In Sec. II we discuss solutions (to arbitrary order in perturbation theory) to the 
density and divergence of the velocity fields in the exact dynamics and non-linear approximations for the Einstein- 
de Sitter universe. Section III is devoted to the statistical description of cosmological fields and the perturbative 
expansion for 1-point cumulants. The diagrammatic approach to perturbation theory is the subject of Sec. IV, where 
we discuss diagrammatic expansions for the variance, skewness, and kurtosis. In Sec. V and Sec. VI we apply the 
formalism developed in the previous sections to the calculation of loop corrections for density and divergence of 
the velocity fields, respectively. The connection of Galilean invariance to the cancellation of infrared divergences is 
considered in Section VII. Section VIII contains our conclusions. Technical material regarding the integrations needed 
to carry out the calculations is given in the Appendices. 
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II. EULERIAN PERTURBATION THEORY 



A. The Equations of Motion 

The N-body problem in which particles interact only through gravity is quite difficult to treat analytically, so some 
simplifications are necessary in order to make progress. The exact evolution equations are equivalent to a hierarchy 
of N coupled integro-differential equations for the n-particle distribution functions in phase space (n = 1,...,N) 
known as the BBGKY hierarchy (see e.g. [6]). By assuming that the particles move without collisions (an excellent 
approximation for weakly interacting dark matter particles) under the influence of the smooth gravitational potential 
of the density field, one can decouple this hierarchy into a mean field equation, the collisionless Boltzmann equation, 
which describes the effective evolution of the 1-particle distribution function. Since we are interested in the study of 
the weakly non-linear regime, we can simplify this further by taking velocity moments of the Boltzmann equation and 
assuming that particle trajectories do not cross (single stream approximation). We also assume the particles have 
negligible thermal motions and move non-relativistically, i.e., that the universe is dominated by cold dark matter 
(pressureless dust). This leaves us with the first two moments (density and velocity fields) as a complete description 
of the system. 

The relevant equations therefore correspond to conservation of mass and momentum and the Poisson equation for 
a self-gravitating perfect fluid with zero pressure in a homogeneous and isotropic universe [6]: 



8t 



+ V-{[l + *(x,r)]v(x,r)} = 0, (2.1a) 



<9v(x, 



h U{t) v(x, r) + [v(x, r) • V]v(x, r) = -V$(x, r), (2.1b) 

u I 

V 2 $(x,r) = ^W 2 (r),5(x,r) . (2.1c) 

Here, x denotes comoving spatial coordinates, the density contrast J(x, r) = p(x, t) / p — 1, with p(r) the mean density 
of matter, v = rfx/rfr represents the velocity field fluctuations about the Hubble flow, % = dlna/dr = Ha is the 
conformal expansion rate, a(r) is the cosmic scale factor, r = J dt/a is the conformal time, $ is the gravitational 
potential due to the density fluctuations, and the density parameter £1 = p/p c = SirGpa 2 /3% 2 . Note that we have 
implicitly assumed the Newtonian approximation to general relativity, valid for non-relativistic matter on scales less 
than the Hubble length aK~ x . We take the velocity field to be irrotational, so it can be completely described by its 
divergence 6 = V • v. By Kelvin's circulation theorem, which follows from taking the curl of Eq. (2.1b), vorticity 
cannot be generated during the non-linear evolution, and any initial vorticity decays as a -2 . This conclusion holds 
true, however, only as long as Eqs. (2.1) are valid; in particular, multi-streaming and shocks can generate vorticity. 
We will refer to Eqs. (2.1) as the "exact dynamics" (ED) (although for the reasons stated above this is a slight abuse 
in terminology), to make a distinction with the modified dynamics introduced by the non-linear approximations to 
be discussed later. Equations (2.1) are valid in an arbitrary homogeneous and isotropic background Universe, which 
evolves according to the Friedmann equations: 

-^-l = --n 2 (T) + -a 2 (r) (2.2a) 

(n-l)n 2 (T) = k-^a 2 (r), (2.2b) 

where A is the cosmological constant, and the spatial curvature constant k = —1,0,1 for Sl tot < 1, Sl tot = 1 and 
&tot > 1 respectively (where Sl tot = SI + Aa 2 /(3W 2 )). 

When the non-linear terms in Eqs. (2.1) are neglected, different Fourier modes evolve independently. Therefore, 
it is natural to Fourier transform the perturbation equations and work in momentum space. Our convention for the 
Fourier transform of a field A(x, r) is: 



A(k, t) = I exp(-ik • x) A(x, r). (2.3) 
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When the non-linear terms in Eqs. (2.1) are taken into account, the equations of motion display the coupling between 
different Fourier modes [5] characteristic of non-linear theories. Taking the divergence of Equation (2.1b) and Fourier 
transforming the resulting equations of motion we get: 



dS(k,r) 



8t 

<9<?(k,r) ~ 3 



0(k,r) = - J d 3 h j d 3 k 2 8 D (k-k 1 -k 2 )a(k,k 1 )6(k 1 ,T)8(k 2 ,T), (2.4a) 



8t 



%(t) 6(k,T) + ^m 2 (T)S(k, T ) = - J d 3 h J d 3 k 2 S D (k-k 1 -k 2 )l3(k,k 1 ,k 2 )6(k 1 ,T)6(k 2 , T ), (2.4b) 



(So denotes the three-dimensional Dirac delta distribution) where the functions 



a(k,k 1 ) = ^, /3(k,k 1; k 2 )= \'<> (2.5) 



k kl riv v i- fc2 ( k i' k 2 
-p-, /?(k,ki,k 2 ) = 22 

IX -y Li f\i -y fx 2 

encode the non-linearity of the evolution (mode coupling) and come from the non-linear terms in the continuity 
equation (2.1a) and the Euler equation (2.1b) respectively. From equations (2.4) we see that the evolution of S(k, r) 
and 6(k, r) is determined by the mode coupling of the fields at all pairs of wave- vectors ki and k 2 whose sum is k, as 
required by translation invariance in a spatially homogeneous Universe [29]. 



B. Perturbation Theory Solutions for Q = 1 and A = 

Equations (2.4) are very difficult to solve in general, since they are coupled integro-differential equations with 
no small dimensionless parameter. The perturbative approach to the problem is to temporarily introduce a small 
parameter A on the RHS of (2.4) to organize a perturbative expansion and then let A — > 1 at the end of the 
calculation [55]. This can be done by taking the fields to be 0(A) and expanding the general solution for S and 6 in 
powers of A. The coefficients of these series can be determined at each order recursively as a function of the lower 
order ones. Ideally one would then like to sum the whole series to recover the solution to the original problem. That 
is of course not possible in practice, so by keeping the first few terms one hopes to get a good approximation to the 
answer. Since the small parameter A is related to the rms density fluctuations, the whole approach is expected to 
break down when these fluctuations become of order one. In this regime, all the terms in the perturbative series are 
of the same order and the perturbation expansion does not make sense without a resummation. 

We consider a matter-dominated Einstein-de Sitter Universe, for which SI = 1 and A = 0. From Eq. (2.2a) we 
obtain a oc r 2 and 3£l% 2 /2 = 6/r 2 . Therefore, Eqs. (2.4) become homogeneous in r and we can formally solve them 
with the following perturbative expansion [5,28,29]: 

oo oo 

£(k,r) = £y(r)<*„(k), £(k,r)=fl(r)$>"(r)Mk), (2-6) 

n=l n=l 

where only the fastest growing mode is taken into account. At small a, the series are dominated by their first terms, 
and since #i(k) = — Si(k) from the continuity equation, <^i(k) completely characterizes the linear fluctuations. The 
equations of motion (2.4) determine S n (k) and n (k) in terms of the linear fluctuations, 

S n {k) = J d 3 qi ... J rf 3 gn fe(k-qi-...-q n )^ s )(qi,...,q n )(5i(qi)...(5i(q n ), (2.7a) 

6 n (k) = - J d 3 gi ...J d 3 q n S D (k - qi - . . . - q n )GW ( qi , . . . , q n )<5! ( qi ) . . . S 1 (q„) , (2.7b) 

where the kernels and Gn ^ are symmetric homogeneous functions of the wave vectors { qi, . . . , q n } with degree 
zero. They are constructed from the fundamental mode coupling functions a(k, ki) and /3(k, ki, k 2 ) according to the 
recursion relations (n > 2, see [5,29] for a derivation): 

F n ((li, . . .,q n ) = ^2 (2^+ 1 3)(w'-T) \^ 2U + ^"O^O-Fn-mCqm+i, • ••><!")+ 2/3(k, ki , k 2 )G n _ m (q m+i , . . .,q n )] , 



(2.8a) 
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n-1 

' (2n + 3)(n ■ 

(2.8b) 



G„(qi, . . .,q„) = ^ (2w+ 1 3)'('w'-T) [ 3a ( k ' kl )- F n-m(qm+i , ■ ■ • , q n ) +2n/3(k,ki,k 2 )G„_ m (q m+ i,...,q. 



(where ki = qi + . . . + q m , k 2 = q m +i + • • • + qn , k = ki + k 2 , and F\ = G\ = 1) and the symmetrization procedure: 

F nH<ll, ■ ■ -,qn) = -7^- Fn ( q7r ( 1 )'---' q7r ( n ))' ( 2 - 9a ) 

7T 

cl s) (qi) • • ->qn) = -7^ Gn ( q7r ( 1 )'---' q7r ( n ))' ( 2 - 9b ) 

7T 

where the sum is taken over all the permutations ir of the set {1, . . . , n). For example, for n = 2 we have: 

pOO/ \ 5 . 1 qi • q2 ( qi , 92 A , 2 (qi ■ g 2 ) 2 

^2 ; (qi>q2) = 7 + 9—— ( — + — j + 7 „2„ 2 . ( 2 - 10a ) 

< z 9i?2 \92 91/ ' 9i9 2 



n {s) ( \ 3 1 qi • q 2 / gi g 2 \ 4 (qi • q 2 ) J ,„,„,, 

G\ qi,q 2 = ■= + TT — + — +- — • 2.10b 

7 2 gig 2 V?2 gi/ 7 

The complexity of the symmetrized kernels increases rapidly with n. For example, the number of terms in F„ 

Is) 

and G n is 134 each for n = 3 and 8523 for n = 4. Explicit expressions for the unsymmetrized kernels i<3 and F4 are 
given in [5]. The perturbation theory kernels have the following properties [5,56] 



1. As k = qi + . . . + q n goes to zero, but the individual q 8 - do not, oc k 2 . This is a consequence of momentum 
conservation in center of mass coordinates. 

(s) Is) 

2. As some of the arguments of F„ or G n get large but the total sum k = qi + . . . + q n stays fixed, the kernels 
vanish in inverse square law. That is, for p ^ qi, we have: 

i#)( qi , . . . , q n _ 2 , p, -p) « Gi s \q u . . . , q n _ 2 , p, -p) (X k 2 /p 2 . (2.11) 



3. If one of the arguments q 8 - of F„ or G„ goes to zero, there is an infrared divergence of the form <\i/q 2 ■ This 
comes from the infrared behavior of the mode coupling functions a(k, ki) and /3(k, ki, k 2 ). There are no infrared 
divergences as partial sums of several wave vectors go to zero. 



C. Perturbation Theory Kernels for Non-Linear Approximations 



When the fluctuations become strongly non-linear, perturbation theory breaks down, and one has to resort to 
N-body simulations to study the subsequent evolution. On the other hand, several non-linear approximation schemes 
have been proposed in the literature which allow analytic calculations beyond the domain of linear perturbation 
theory. These are approximations in the sense that they replace the Poisson equation by a given ansatz which is true 
only in linear theory for the exact dynamics [52], and therefore they are neither exact nor asymptotic to the exact 
solution beyond linear order (except for special cases). In this section we derive the perturbation theory kernels in 
Fourier space for the Einstein-de Sitter model for these approximations. Similar analysis has been done in real space 
for second order perturbation theory in [52] and up to third order in [51]. 
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1. Zel'dovich Approximation (ZA) 



In this approximation [46,57], the motion of each particle is given by its initial Lagrangian displacement: the 
dynamics of fluid elements is therefore governed purely by "inertia". In Eulerian space, this is equivalent to replacing 
the Poisson equation by the ansatz [52]: 



3H(t) 



-V$(x,t) 



(2.12) 



which is the relation between velocity and gravitational potential valid in linear theory. The important point about 
ZA is that a small perturbation in Lagrangian fluid element paths carries a large amount of non-linear information 
about the the corresponding Eulerian quantities, since the Lagrangian picture is intrinsically non-linear in the density 
field. This leads to non-zero Eulerian perturbation theory kernels at every order. ZA works reasonably well as long 
as streamlines of flows do not cross each other. However, multistreaming develops at the location of pancakes leading 
to the breakdown of ZA [57]. The equations of motion in Fourier space are: 



dS{k,T) 
dr 



(k,r) 



h / rf 3 fc 2 fe(k-k! -k 2 )a(k,ki)£(ki,r)£(k 2 ,r) 



(2.13a) 



S0(k,r) U(t) 



h / d d A 2 <y i? (k-ki-k 2 )^(k,ki,k 2 )(9(ki,r)(9(k 2 ,r) 



(2.13b) 



These equations, together with the perturbative expansion (2.6), lead to the recursion relations (n > 2) 

n-l 



^ n z (qi> • • - ,qn) = G m(qi> • ••><!« 



»( k ,ki) F Z / 

r n-m\'i™ + l 



/3(k,ki,k 2 ) z 

n — ml1 m +l' * * * ' q^ 



n(n — 1) 



(2.14a) 



r <z t \ _\r K r iz i N^(k, ki,k 2 ) z 

(.qi, • • • , q«j — >. (j- m (.qi, • • • , qmj — -, 77 — cr„_ m (q m+ i, . . . , q„j, 

z — ' (n — 1 



(2.14b) 



In this case, we can factor out the wave vector dependence in the denominator by defining: 



F n (qi, • • - ,qn) 



— 2 -j^f (qi, . . .,q n ) 

1\ ■ ■ -In 



(2.15a) 



G>f(qi, . • - ,q n ) 



91 • • - 9 

where, again, k = qi + . . . q n . The recursion relations become 



2 „2 $n (qi, • • • , q«) 



(2.15b) 



J-^n (qi, • • - ,qn) = ^ ^m(qi) • • - ,q m )(k -ki)^f_ m (q m+ i, . . .,q n ) + & 2 £,f(qi, • . .,q„ 



(2.16a) 



2 («- i)£n (qi, • • -,qn) = X] ^m(qi> • • -,q m )(ki -k 2 )£;f_ m (q m+ i,...,q„), (2.16b) 

m = l 

where t/f = 1 and ^"^(q) = q 2 ■ From these equations it is clear that all the wave- vector dependence in the denominator 
of the kernels has been factored out in (2.15) and is not present in Q or T . In fact, after symmetrization, the density 
field kernel takes the simple form [50]: 

pz(*)r„ 1 (k- qi ) (k • q n ) 

t n ^(qi, . . .,q n ) - — — 2 — ... — - 2 . (2.17) 

9l In 
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Note that this expression means that in ZA, property (1) for given in the previous section is changed to oc k n . 

So, in ZA, higher order kernels are underestimated for large scales compared to ED. For one-dimensional perturba- 
tions, the ZA becomes exact, and therefore its symmetrized kernels (where dot products are replaced by ordinary 
multiplication) agree with the symmetrized kernels of the exact dynamics given in Section II B. In 3 dimensions, for 
n = 2 we have (compare (2.10)): 

p z(<), \ i . i qi • q2 f qi . 92 A , 1 (qi • q.2) 2 ,„ 18 , 

F 2 (qi,q2 = 77 + 77 1 + o — — > 2.18a 

2 2 q x q 2 \q 2 qi J 2 q(q 2 

-< z (. s )t„ „ \ _ 1 qi • q 2 fn | iA | 1 (qi • q2 ) 

2 qiq 2 V92 qi J 2 q\q\ 



G 2 V ; (qi,q 2 ) = ^— — — + — ) + 7;— -T7T- ■ (2.18b) 



In this case, the number of terms in Fn^ and increases more slowly than in the exact dynamics kernels. There 

are 26 (16) terms in F^ (s) (Gg (s) ) and 237 (127) in F% (s) (G^ (s) ). 

2. Linear Potential Approximation (LPA) 

In this approximation [48,49], the gravitational potential is substituted by its linear value. Therefore, the full 
Poisson equation (2.1c) is replaced by: 

V 2 $(x,t) = ^ 2 (t)Mx,t), (2.19) 

where Si(x,t) = a(r)Ji(x) is the solution to the linearized equations of motion. The idea behind this approximation 
is that since $ oc S/k 2 , the gravitational potential is dominated by long- wavelength modes more than the density 
field, and therefore ought to obey linear perturbation theory longer. In fact, N-body simulations show that $ evolves 
much more slowly than the density field [48]. The equations of motion in Fourier space are: 



dS(k,r) 



8t 



k,r) = - J d 3 h J rf 3 fc 2 fe(k-k! -k 2 )a(k,ki)£(ki,r)£(k 2 ,r), (2.20a) 



U(t) e(k,T) + ^n 2 (r)S 1 (k,T) = - J d 3 h J rf 3 fc 2 fe(k-k! -k 2 )/3(k,k 1 ,k 2 )^(k 1 ,r)^(k 2 ,r), (2.20b) 



8t 

leading to the following recursion relations (n > 2) 

n-l 

i^ p (qi,...,q„) = J2 G™ P (qi,...,q 



r ^Mi)piP / n , , 2/3(k,k!,k 2 ) ^ LP 

n n(2n + 1) 



(2.21a) 



G^ P (qi, • • -,qn) — V] G L f(qi, . . .,q m ) ^ ' 1; ' 2 ^ G^ p m (q m+ i, . . .,q n )- (2.21b) 
' (2n + 1) 

m = l v ' 

Making the analogous definitions to (2.15), we find: 

n-l 

n^ P (qi, . . . , q„) = ^ ^ P (qi, • • • , q m )(k • k^;^ (q m+1 , . . . , q„) + A 2 # P (qi, . . . , q„), (2.22a) 

m = l 

n-l 

(2n + l)^ P (qi,...,qn) = V Q L J (qi , . . . , q m ) (ki • k 2 )Q L n P _ m (q m+1 , . . . , q„) . (2.22b) 
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Note that (2.22a) is exactly equivalent to its Zel'dovich approximation counterpart (2.15a), whereas (2.22b) has a 
different numerical coefficient from (2.16b). For n = 2 we have: 

P tP( s )/ x 1 . 7 qi • q 2 (qi q 2 \ , 1 (qi • q.2) 2 ,„ n 
F 2 qi , = 77 + 777T h— +- — — , 2.23a 

2 20 q x q 2 \q 2 qi J 5 q(qi 

n LP( s ), , l qi • q2 f qi .92 V 2 (qi ■ g 2 ) 2 ,„„„,> 

G 2 ; (qi,q 2 ) = t— — [ — + — +7 „2„2 • ( 2 - 23b ) 
9l?2 \92 91/ q 1 q 2 

The number of terms in the symmetrized kernels for n > 2 is equal to the ZA case. 

3. Frozen Flow Approximation (FFA) 

In this approximation [47] the velocity field is assumed to remain linear. Therefore, we have the relation: 

0(x, r) = h (x, r) = -ft (r)*i (x, r). (2.24) 

Streamlines are kept frozen to their initial configuration, therefore the dynamics remains forever in the single-stream 
regime. The equation of motion for S(h, r) is: 



dS{k,T) 



n(T)S 1 (k,T)=n(T) J d 3 h J d 3 A 2 <fo(k-ki -k 2 )a(k,ki)£i(ki,r)£(k 2 ,r), (2.25) 



8t 

where Ji(k, r) = a(r)Ji(k). This leads to the perturbation theory kernels (n > 2): 

p-FF / n x _ 1 (k-qi) [(k-qi) • q 2 ] (q n • q n ) 

(qi, • • • , qn) - —f ~2 ~ 2 ••• ~ 2 , (2.26a) 

1\ 1 2 In 

G^ F ( qi ,...,q n )=0. (2.26b) 

For n = 2, we have: 

^ F(s) (q 1 ,q2) = | + i^f^ + ^. (2.27) 
2 4 gig 2 \q 2 q x J 

In this case, the number of terms in F„ F ^ is 16 and 125 for n = 3 and n = 4, respectively. 

It is worth noticing that all these non-linear approximations involve neglecting "a coupling" (for ZA and LPA 
and also "/? coupling" for FFA) in the velocity divergence kernels, which comes from the fact that none of these 
approximations solves Poisson equation. This leads to the result that both kernels in these approximations do not 
involve scalar products of wave- vectors in their denominators. This fact will significantly simplify the calculation of 
1-point cumulants. 

III. STATISTICAL DESCRIPTION OF FLUCTUATIONS 

The starting point for a statistical description of fluctuations in cosmology is the "Fair Sample Hypothesis" [6,58]. 
This asserts that fluctuations can be described by homogeneous and isotropic random fields (so that our Universe 
is a random realization from a statistical ensemble) and that within the accessible part of the Universe there are 
many independent samples that can be considered to approximate a statistical ensemble, so that spatial averages are 
equivalent to ensemble averages ("ergodicity"). 

A full description of a set of random fields is provided by their joint probability distribution functional. However, 
except for special cases, the determination of this functional cannot be carried out exactly [20], and in practice one 
considers only the simplest statistical tools which describe some particular statistical property of the fluctuations. A 
very important set of functions are the moments of the probability functional. The (equal-time) p th -order, p-point 
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moment for the density field is defined by the ensemble average (denoted by angle brackets) of the product of p fields 
at p different points (in the following, we give formulas in real space, but analogous expressions hold in Fourier space; 
see, e.g., [59] for a general discussion): 



M x i> 



(S(xi)...S(x p )), 



(3.1) 



where jj, p (xi, . . . , x p ) = £(xi, . . . , x p ) is usually referred to as the p-point correlation function. However, in general 



, x p ) with p < 



one may have p-point moments of order r (p, r (xi 

fields coincide with each other. Note that in (3.1), evaluation of all the fields at the same instant is understood 
moments can be obtained from the moment generating functional Z( J) by functional differentiation: 



when some of the points of evaluation of the 

The 



-iyspz[j{x)] 



SJ(xi) . . . SJ(x p ) 



j=o 



where: 



Z[J{x)] = <^exp [- / d 3 xS{x)J{x)]^ = / VS{x)P[S{x)]exp[- j d 3 xS{x)J{x)] 



(3.2) 



(3.3) 



and P[J(x)] is the one-point probability distribution functional of the density field fluctuations. 

Under very general conditions, when all the moments exist, they characterize the underlying probability functional 
[60,61]. However, they are not the only set of functions with this property, or even the best. One inconvenience is 
that each moment contains correlations already present at lower orders. For that reason, one defines another set of 
functions, the cumulants K p , which specify the non-trivial correlations at each order. These are given in terms of the 
recursion formula (together with the condition \i\ = K\ = (8) = 0): 



M x i> 



Kp (Xl , 



E 

ul a =l 



(3.4) 



where I = {xi, . . . , x p }, \I a \ is the number of elements of the subset I a , and we sum over all the possible partitions 
of the set of arguments excluding the identity partition. Since all the trivial correlations have been taken away, 
the cumulants vanish when any subset of X is removed to infinite separation. The p th order, p-point cumulant 
K p (xi, . . . , x p ) is usually referred to as the connected correlation function (in analogy with the connected Green's 
functions in quantum field theory [62]) and also denoted as £ c (xi, . . . , x p ) or (S(xi) . . . S(x p )) c . For isotropic and 
homogeneous random fields, the cumulants depend only on the relative distances |x 8 - — Xj \ for i, j = 1, . . . ,p (i ^ j). 
In Fourier space, homogeneity implies that the configuration in k-space is closed, that is Yli=i kj = 0, whereas isotropy 
implies that the cumulants are invariant under reorientation of this closed configuration. For a Gaussian probability 
distribution functional, K p = for p > 2. 

The cumulant expansion theorem states that the cumulant generating function is simply lnZ[J(x)] (the proof is 
easily done by induction, see e.g., [63]), therefore: 



K p (Xl , 



-1)P6P In Z[J(x)} 



_ SJ(xi) . . . SJ(x p 



(3.5) 



j=o 



which gives a practical way of calculating cumulants when the probability functional is known. 

In this paper, we adopt the conventional assumption that the primordial fluctuations are Gaussian distributed. 
While not universal, this is a generic prediction of the simplest inflationary models, in which Si is linearly related 
to an essentially free scalar field (the fluctuation of the inflaton field) in its vacuum state, described by a Gaussian 
wavefunctional. For Gaussian fluctuations, lnZ[J(x)] is quadratic in J(x), and therefore K p vanishes for p > 2 for 
the linear fluctuations. Since the equations of motion of a self-gravitating fluid are non-linear, higher order (p > 2) 
cumulants will be induced by dynamical evolution, driving the probability functional away from Gaussianity [6]. 
Once this occurs, the cumulant generating functional is no longer known a prion, and (3.5) cannot be employed 
to evaluate the cumulants (see, however, [19] for a calculation of the cumulant generating functional when the rms 
density fluctuation is vanishingly small). 

In this work we focus on the 1-point cumulants of cosmological fields induced by non-linear gravitational dynamics. 
The linear fluctuations are characterized by the second order cumulant, the variance af = {Sf) c , or equivalently, by 
the linear power spectrum Pi(k, r), defined by: 



Mk)Mk'; 



S D {k + k')Pi{k,r) 



(3.6) 
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From (3.4), for Gaussian initial conditions we have: 

^i(ki) . . .Ji(k 2p )^ = ^- £ (M k 7T(i))M k 7r(2))) c • • • (^i(k 7r (2p-i))<5i(k 7r(2p) )^ (3.7a) 
frfi(ki) .. .Ji(k 2p+1 )) = (3.7b) 



Due to homogeneity, one-point cumulants are independent of the spatial coordinate, and in Fourier space we have: 



(<F(x)) c = {SP) C = J d d h . . .d'kp^SQm) . ..S(k p )) c . (3.8) 

Thus the p th order, one-point cumulants are integrals of the p-point spectra over the p wave- vectors. According to 
the perturbation expansion (2.6), (S p ) c can be written as: 

OO p 

r=p-l C(2r,p) i = l 

where we sum over the arrangements C(2r,p) of p positive integers that sum up to 2r = 2p — 2,2p, . . . (because of 
the assumption of gaussian initial conditions, 2p — 2 is the minimum order in Si that we need to get a connected 
contribution [7]) and Ci denotes the i th component of a particular arrangement (i.e. an integer between 1 and 2r— p+l). 
In Eq. (3.9), terms labeled by r are contributions to (6 P ) C of order 6f r . For example, for p = 2,3,4 and working to 
0{6\ ) we have (putting contributions of same order within square brackets; here and below superscript (n) denotes 
the n-loop contribution): 

2 = /A 2 \ -,2/i , , (2) 4 , (3) 6\ , n( 10 



<r z = (6 Z ) C = ai (l + s^al + s^a\ + s^J + 0{a^), (3.10a) 

S W = [2(S 1 S 3 ) c + (S? 2 ) c ] , (3.10b) 

S < 2 ) = [2(S 1 S 5 ) c + 2(S 2 S 4 ) c + (61) c ] , (3.10c) 

= [2(S 1 S r ) c + 2(S 2 S 6 ) c + 2(S 3 S 5 ) c + (Sl) e ] , (3.10d) 

{S 3 ) c = 3(6fa) e + [3(SlS 4 ) e + 6(W 3 ) C + (<*f )J + [HS'fS 6 ) c + 6(Ws) c + 6(6 1 6 3 6 4 ) c + 3(6*6 4 ) e 

+oJs 2 Sl)]+O(al ), (3.11) 



(6% = [4(6*63) e + Q(S!6l) c ] + [4(6f6 5 ) c + I2(6f6 2 6 4 ) c + 6(6*6%) e + I2(6 1 6l6 3 ) c + (6 4 2 ) c ] + 0{<r\% (3.12) 

where the numerical factors come from the multinomial expansion of 6 P in 6 n 's, and at is the linear rms density 
fluctuation given in terms of the linear power spectrum by: 

aj = (6l) c = J dPkP^k^) (3.13) 

In view of (3.9), a convenient way to measure the departure from gaussianity due to dynamical evolution is to 
introduce the S p parameters (p > 3) defined by rescaling the 1-point cumulants [5]: 



(6 2 



Sp = -^LJ^. (3.14) 



Therefore, the S p parameters are constants and independent of the normalization of the linear power spectrum to 
lowest order in perturbation theory. Note that since the S p 's are constructed directly from one-point cumulants, each 
of these parameters contains statistical information independent of the others. S3 is a measure of the skewness of 
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the probability distribution (a positive value generally indicates that the upper tail of the distribution is the higher), 
whereas S4 is a measure of the kurtosis (i.e., how "flat" or "peaked" the distribution is compared to a gaussian) [61]. 
From (3.10) we have 

S3 = 4 0) + 4 1) ^ 2 + S { i ] 4 + O(af), (3.15a) 
S 4 = 5l 0) + st ) aj + 0{a A l ), (3.15b) 

where: 

5 (o) 3(ffi 2 ) e (i) 3(ffi 4 ) e + 6(<W 3 ) e + (S%) c _ 2 (i) 5 (o) f3 16a) 

3 4 ' 3 6 3 ' V * / 

5 ( 2 ) = 3(^4) e + 6( W 5 ) e + 6( W 4 ) e + 3(^ 4 ) e + 3(^ 3 2 ) e _ 2g(1)5 ( 1} _ , (1))2 + 2s ( 2) \ 5 (o ); (3 16b) 

a i ^ ' 

c(Q) _ 4(ffi 3 ) e + 6(<^ 2 2 ) e ^(1) _ 4(ffrf 5 ) e + 12(^ 4 ) e + 6(ffig) e + 12(<M 2 2 ^)c + (^|)c 3j(1)5 (o) (g lgc) 

and s^ 1 ) and s 1 - 2 - 1 come from the next-to-leading order corrections to {S 2 ) c in the denominator of Eq. (3.14) (see 
Eq. (3.10a)). To calculate the necessary ensemble averages, we use (2.7) to relate S n to the linear fluctuations 
and then ergodicity to relate averages over the present probability distribution to averages over the initial gaussian 
ensemble. In general we have: 

00 „ „ 

( SP )c= E E M 3 ?i--7 rf M^V---^j---^;W (3.i7) 

r=p-l C(2r,p) c 

The next step is to reduce the moments of linear fluctuations in (3.17) according to their gaussian nature given 
by (3.7) and then take into account only the connected contributions to the integrals. The complexity of this procedure 
can be easily handled by using diagrammatic techniques. 

IV. DIAGRAMMATIC APPROACH TO PERTURBATION THEORY 

A. Diagrammatic Rules 

A systematic framework for calculating correlations of cosmological fields in perturbation theory has been formulated 
using diagrammatic techniques similar to those in quantum field theory [62] and statistical mechanics [63]. In this 
approach [5,56], contributions to p-point cumulants of the density field come from connected diagrams with p external 
(solid) lines and r = p — l,p, . . . internal (dashed) lines. The perturbation expansion (3.17) leads to a collection of 
diagrams at each order, the leading order being tree-diagrams, the next to leading order 1-loop diagrams and so on. 
On the other hand, Eqns. (3.10), (3.15) indicate that the perturbation series for the one-point cumulants can be 
viewed as an expansion in powers of cr|. Thus, in analogy with quantum field theory, cr| is the effective 'coupling 
constant' for the loop expansion, and we expect perturbative results to be sensible for small coupling, af 1. 

In each diagram, external lines represent the spectral components of the fields we are interested in (e.g., S(h, r)). 
Each internal line is labeled by a wave- vector that is integrated over, and represents a linear power spectrum P\{q, r). 
Vertices of order n (i.e., where n internal lines join) represent a n th order perturbative solution S n , and momentum 
conservation is imposed at each vertex. Figure 1 shows the factors associated with vertices and internal lines. To find 
the contribution of order 2r to the p-pomt spectrum of the density field proceed as follows: 

• Draw all distinct connected diagrams containing p vertices (with external lines labeled by ki . . . k p ) joined by r 
internal lines. Two diagrams are distinct if they cannot be deformed into each other by moving the vertices and 
lines without cutting any internal lines (sliding lines over other lines is allowed in the rearrangement process). 
For each of these diagrams: 
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1. Assign a factor of (k 8 - — qi — . . . — q n )-Fn (qi , . . . , q n ) to each vertex of order n and external momentum k 8 - 
(i = 1, . . .,£>). For the arguments of , we use the convention of assigning a positive sign to wave- vectors 
outgoing from the vertex. 

2. Assign a factor of P\{qj, t) to each internal line labeled by q^. 

3. Integrate over all q^ (j = 1, . . . , r). 

4. Multiply by the symmetry factor of the graph, which is the number of permutations of linear fluctuations 
(Si's ) that leaves the graph invariant. 

5. Sum over distinct labelings of external lines, thus generating p\/(ni\ . . .ri2 r - p +i!) diagrams (where n 8 - 
denotes the number of vertices of order i) . 

• Add up the resulting expressions for all these diagrams. 

To calculate 1-point cumulants, we have to integrate further over the k 8 (i = 1, . . .p) according to (3.8). These 
integrations are trivial because of the presence of delta functions given by rule 1, so we are left only with integrations 
over the qj's. Also, since the p\/(ni \ . . .ri2 r - p +i!) diagrams generated by rule 5 become equal contributions when the 
integration over external lines is performed, to find the contribution of order 2r to the p th -order 1-point cumulant of 
the density field we replace rule 5 by the following: 

• 5a. Integrate over k 8 (i = 1, . . .p) and multiply by the multinomial weight p\/(n-\ \ . . .ri2 r - p +i!). 



S D (k - qi - . . . - q n )-Fn (qi, • • • , q n ) 
Pi(q,r) 

FIG. 1. Diagrammatic rules for vertices and internal lines for density field fluctuations. Equivalent rules hold for the velocity 
divergence replacing F„ 3 ^ by . 




A few comments about these rules are in order. Rule 1 basically assigns a perturbation theory kernel of order 
n to each vertex of that order. Rule 2 is a consequence of the assumed gaussian nature of the linear fluctuations, 
and describes the fact that the Si 's coming out from each vertex "interact" with each other only through two-point 
functions. Rule 3 (together with rule 1) says that contributions to each external line with wave- vector k 8 come from 
the mode coupling of the linear fields at wave vectors q^ whose sum is k 8 , as required by translation invariance. Rule 4 
counts the number of ways we can join Si 's into two-point correlations. Rule 5a comes from the multinomial expansion 
of S p into 8 n 's according to the perturbation expansion (2.6). The same diagrammatic rules hold for the divergence of 

Is) Is) 

the velocity field, replacing F„ by G n . Using both sets of rules, one can construct diagrammatic expressions which 
represent contributions to arbitrary cross-correlations between density and divergence of the velocity fields. Finally, 
we note that these rules are somewhat unusual from the viewpoint of field theory, since new vertices appear at each 
order of perturbation theory. 

It is convenient to introduce some bookkeeping notation in order to be able to refer to different diagrams. 
We shall denote by X> 8 - j the amplitude given by the above rules for a diagram representing the contribution to 
(Si(qi) . . .Sj(q p )) c (e.g., see Fig. 2), and define 

^ _ Jd<Tj((li)...Sd<Tj(q p )Vi... j 

fdoj( qi )...fdoj( qp ) ' (4 ' ij 

where the bar denotes the average of diagram amplitudes over the initial rms fluctuations, and: 

daj( q ) = d 3 qPi( q ,T). (4.2) 

We now turn to the explicit construction of the diagrammatic series for each cumulant to be considered. 
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B. Loop Expansion for a 2 



The diagrams corresponding to the perturbative calculation of the variance (2nd order 1-point cumulant) are given 
in Figs. 2, 3, 4 and 5. According to the above rules, their contribution is: 



-r 2 _ 2 

U = (Tii 



where: 



J ^'(qi) J d<rj{<l2) [V22 + 2>3i] + J ^'(qi) - J do-jiqs) [v£ 3 + Vg + V 42 + 2? 5 i] 
J ^(qi) ...J d<7 2 (q 4 ) [v{ 4 + V{{ + V{ 3 + Vg + V 62 + V 



2? 22 = 2^ 2 (s) ( qi ,q 2 )^ 2 (s) (-qi>-q2) = 2[^ 2 (s) ( qi ,q 2 )] 2 , 

2>31 = 6F3 (,) (qi,-qi,q 2 ), 



(4.3) 



(4.4a) 
(4.4b) 



are the one-loop contributions (see Fig. 2), 

^3 = 6^3 (s) ( qi ,q 2 ,q3)^3 (s) (-qi>-q2,-q3) 



9^ 3 (S) (qi , -qi , q2)^ 3 (s) (-q2 , q 3 , -q 3 ) 



(4.5a) 
(4.5b) 



V 42 = 24F 4 (,) (qi , - qi , q 2 , q3)F 2 (,) (-q 2 , -q 3 ) 



V 51 = 30-F 5 w (qi,-qi,q 2 ,-q 2 ,q 3 ), 
correspond to the two-loop diagrams (see Fig. 3), and 

£>4 4 = 24_F 4 (s) (qi,q 2 ,q3,q 4 )_F 4 (s) (-qi,-q 2 ,-q3,-q 4 ) 
£>44 = 72_F 4 (s) (qi,-qi,q 2 ,q 3 )_F 4 (s) (-q 2 ,-q3,q4,-q 4 ) 



(4.6) 
(4.7) 

(4.8a) 
(4.8b) 



V 1 



V 



II 



120_F 5 (s) (qi, -qi, q 2 , q 3 , (li)F { 3 s) {-q_ 2 , -q_ 3 , -q 4 ) 
90_F 5 (s) (qi, -qi, q 2 , -q 2 , q 3 )-F 3 (s) (-q3, q4, -q4), 



(4.9a) 
(4.9b) 



X»62 = 180_F 6 (s) (qi,-qi,q 2 ,-q 2 ,q 3 ,q 4 )-F 2 (s) (-q3,-q4) 



(4.10) 



71 



210Ff ; (qi, -qi, q 2 , -q 2 , q 3 , -qs, q4) 



(4.11) 



describe the 3-loop corrections (see Figs. 4 and 5). 

The V22 diagram has a symmetry factor of 2, corresponding to the two ways in which the 4 Si's (2 in each vertex) 
can be paired. The T>si diagram has a symmetry factor of 3, since there are 3 ways in which we can choose one Si 
out of 3 Si's in the third order vertex to join with the only Si coming out from the first order vertex. There is also 
a multinomial factor of 2 in this case. The numerical factors in the 2-loop and 3-loop expressions are derived in a 
similar fashion. 

Using Eq. (4.1), we can write (see Eq. (3.10)): 



.(2) 



.(3) 



22 

w 



■2>31, 

-vU + v 4 



v{ 4 + v{{ + vi 3 



■v 51 , 



■T>7 



(4.12) 
(4.13) 
(4.14) 



which express loop correction coefficients in terms of averages of diagram amplitudes over the linear rms density 
fluctuations. 
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qi 



qi 



q2 



q2 



(2>n) (V 22 ) (D31) 

FIG. 2. Diagrams for a 2 up to one loop. See Eqs. (4.4) for diagrams amplitudes. 
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q2 
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qs 
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qi 
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q2 



qs 
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qi 

'*~* q2 



qs 



qi 

q2 n 
*- i 

t > ' qs 



(2? 4 



(2>5 



FIG. 3. Two-loop diagrams for <r 2 . See Eqs. (4.5), (4.6) and (4.7) for diagrams amplitudes. 



qi 
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q 2 ^ 
- ~ - \ 

q3 
q~4 



qi 


q4 


q 2 




( t"*- 








~~ - ► - 




qs 





qi 
*- 
q2 
»- 



^ | qa 



q4 



(^ 4 ) (2>62) 
FIG. 4. Three-loop diagrams for <r 2 corresponding to D44 and De2- See Eqs. (4.8) and (4.10) for diagrams amplitudes. 



qi 



s q 2 

■ i 

q4 



qi 
, *- 

I 2 ) 



q4 

( i 



qi 
, »- ^ 

/ q2 \ 

,' ' *■ - \ 
1 ' q3 n i 
>i ». 1 

\\ ( i /, 



q 4 



(^3) (Vl£) (V 71 ) 

FIG. 5. Three-loop diagrams for a 2 corresponding to T>$z and X>7i- See Eqs. (4.9) and (4.11) for diagrams amplitudes. 



C. Loop Expansion for S3 

From the diagrammatic rules and Fig. 6 we have: 

^211 = 6 J F 2 (s) (q 1 ,q 2 ), 



(4- 



for the tree-level contribution. 
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q 2 ^__ 

(2>21l) 

FIG. 6. Tree-level diagram for (5 3 } c - The corresponding amplitude is given by Eq. (4.15). 



qs T + — '^2 + 

(^222) (2>£ 21 ) (2>&) (^11) 



FIG. 7. One-loop diagrams for (S 3 ) c - The corresponding amplitudes are given in Eq. (4.16). 

For the 1-loop diagrams in Fig. 7 we get: 

V 222 = 8F 2 is \ qi , q2 )F 2 is \ q2 , q3 )F 2 is) (q 3 ,- qi ), (4.16a) 

2>32l = 36^ 3 (s) (qi,q 2 ,q 3 )i ? 2 (s) (qi,q2), (4.16b) 

VH, = 36F 2 is \ qi , q2 )F^ s \- q2 , q3 ,- q3 ), (4.16c) 

2>411 = 36 J F 4 (s) (qi,-qi,q 2 ,q 3 ). (4.16d) 

From Figs. 8, 9, and 10 we get the 2-loop contribution to the third order cumulant: 

Vi 32 = 108^3 (s) (q 1 ,q 2 ,q 3 )^3 (s) (-qi,-q2,q4)i ? 2 (s) (q3,q4), (4.17a) 

Vg 2 = 108^ 3 (s) (q 1 ,q 2 ,q 3 )^3 (S) (-q3,q4,-q4)i ? 2 (s) (qi,q 2 ), (4.17b) 

^332 = 54^ 2 (s) (q 1 ,q 2 )^3 (s) (-q 1 ,q 3 ,-q3)i ? 3 (S) (-q2,q4,-q4), (4.17c) 

V{ 22 = 72^ 4 (s) (q 1 ,q 2 ,q3,q 4 )i ? 2 (s) (qi,q 2 )i ? 2 (s) (q3,q4), (4.18a) 

V{[ 2 = 144^ 4 (s) (q 1 ,-q 1 ,q 2 ,q 3 )^ 2 (s) (q 2 ,q4)i ? 2 (s) (-q3,q4), (4.18b) 

V{ 31 = 144 J F 4 (s) (q 1 ,q 2 ,q3,q4) J F 3 (S) (-qi>-q2,-q3), (4.19a) 

2>& = 2i6^ 4 (s) (-qi>-q2,q4,-q4)i ? 3 (s) (qi>q2,q3), (4.i9b) 

ViH = 216^ 4 (s) (q 1; -q 1; q 2 , q 3 )^ 3 (S) (-q3, q 4 , -q 4 ), (4.19c) 

V{ 21 = 360F 5 (,) (qi, -qi, q 2 , q 3 , ^F^ (q 2 , q 3 ) , (4.20a) 

2>52l = 180i ? 5 (s) (-q2,q3,-q3,q 4 ,-q 4 )i ? 2 (s) (qi,q 2 ), (4.20b) 

V 611 = 270 J F 6 (s) (q 1 , - qi , q 2 , -q 2 , q 3 , q 4 ). (4.21) 
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(^332) (^332) (^332) 

FIG. 8. Two-loop diagrams for (S 3 ) c corresponding to X>332. See Eqs. (4.17) for diagrams amplitudes. 



qi 



,<q2 

■*<13 



q4 



(v{ 2 



qi 

✓ q2 
1 .+ 



iq4 



qs 



q2 



♦r^-'^qs 

qi 



qi 



q2 



(^21) 



qs 
q4 



FIG. 9. Two-loop diagrams for (5 3 } c corresponding to D422 and D521. See Eqs. (4.18) and (4.20) for diagrams amplitudes. 



/, -q 2 .V q3 



q 4 ^. 



qi 



q 4 

( ) 



s q2 
q3^„ 



qi 

, *• s q2 . 



q 4 



q3 ^-^v 



qi 

q 2 ^ qs . 



q4~ 



(^31) 



(v 



611, 



(*>43l) (^431, 
FIG. 10. Two-loop diagrams for (S 3 ) c corresponding to D431 and Den- See Eqs. (4.19) and (4.21) for diagrams amplitudes. 

Therefore: 



<* 3 



J daj^)... J da 2 i ( (l3 )[v 222 +Vi 21 +Vi I 21 +V 411 ] + J do-fa 



J daj(<l4) ^332 + V 332 + D332 + X>4 2 2 + P422 + T>{ 31 + V{ 31 + + T>\ 21 + + I>611 



+0(<ri 



(4.22) 



which gives the contributions to the 3rd-order cumulant up to 2-loop corrections. Here, the tree-level diagram X>2ii 
(see Fig. 6) has a symmetry factor of 2! due to the permutations in the second order kernel. The one-loop diagrams 
(see Fig. 7) have symmetry factors of 8 for X>222 (2! for each second order kernel), 6 for X>3 21 (3 to choose one S\ to 
connect the third order with the first order kernel, 2! from the remaining permutations with the second order one), 
6 for P321 (3 from choosing one S\ in the third order kernel, 2! from the second order one) and 12 for X>4n (6 from 
choosing two Ji's out of four in the fourth order kernel and 2 from the remaining permutation). Similar analysis gives 
the numerical coefficients in the two-loop diagrams. 
We can now write Eq. (3.16a) as: 



2>211 
2? 225 



V 1 

^321 



V II 

^321 



■ X>4n 



(4.23a) 
(4.23b) 
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and Eq. (3.16b) as: 

S3 ' = T>1 32 + + ^332 + ^422 + ^422 + ^431 + ^431 + ^431 + ^521 + ^521 + ^611 — 2s 1 " - 1 S3 ' 

_(( s (l))2 + 2s (2)^(0)^ (4 24) 



D. Loop Expansion for S4 

According to the diagrammatic rules in the case of the fourth order one-point cumulant we have: 



2>22ii = 48F 2 W (q 1 ,q 3 )F 2 W (-q 3 ,q 2 ), (4.25a) 

2?3111 = 24 J F 3 (s) (q 1 ,q 2 ,q3), (4.25b) 

for tree-level diagrams (see Fig. 11). 

For one-loop diagrams (see Figs. 12, 13 and 14): 

2? 2222 = 48^ s) (q 1 ,q 2 )^ 2 (s) (-q 2 ,q 3 )^ 2 (s) (q3,q4)i ? 2 (s) (q4,-qi), (4.26) 

V\ T11 = 288^ 2 (s) (q 1 ,q 2 )^ 3 (s) (-q2,q3,q4)i ? 2 (s) (qi,q3), (4.27a) 

Vg 2l = 288^3 (s) ( qi ,q 2 ,q3)^ 2 (s) (qi,q 2 )^ 2 (s) (-q3,q4), (4.27b) 

V 1 ^ = 288^ 2 (s) (q 1 ,q 2 )^ 3 (s) (-qi>q4,-q4)i ? 2 (s) (-q2,q3), (4.27c) 

Vi 311 = 216^ 3 (s) ( qi , q 2 , q 3 )^ 3 (s) (q4, -q 2 , -q 3 ), (4.28a) 

Vg n = 216^ 3 (s) ( qi , q 2 , q 3 )^ 3 (S) (-qi> q 4 , -q 4 ), (4.28b) 

V{ 211 = 288 J F 4 (s) (q 1 ,q 2 ,q3,q 4 ) J F 2 (s) (qi,q 2 ), (4.29a) 

VH U = 576^ 4 (s) ( qi ,- qi ,q 2 ,q3)^ 2 (s) (-q3,q4), (4.29b) 

2> 5m = 240 J F 5 (s) (q 1 , - qi , q 2 , qa, q 4 ), (4.30) 



qi qi 



' \q2 

q3T + q3T \ 



q2 



(2>221l) (2>311l) 
FIG. 11. Tree-level diagrams for (S 4 ) c - The corresponding amplitudes are given in Eq. (4.25). 
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Therefore: 



{S 4 ) c = J da'Ufn) . . . J d<7|(q3) [d 22 h + 2)311i] + J d<7|(qi)...y d<7|(q 4 ) [d 2222 + D^ 221 + 2? ; 



) /J 

y 3221 



+^3221 + ^3311 + ^3311 + ^4211 + ^4211 + ^5 



5111 



0(a 



10\ 



so we can write Eq. (3.16c) as: 

AO) 



SI 

s: 



V, 



+ V 31 



(i) 



V 



2222 + 2^3221 + ^3221 + ^3221 + ^3311 + ^3311 + ^4211 + ^4211 + ^5 



5111 



(4.31) 



(4.32a) 
(4.32b) 



qi qi 



. — / _ N - — , 

i i / \ i 

i i / \ i 

q2f rq4 + q3» » q2 + q3f 

i i \ / i 

'< qs ! \ / q4 





q4 






qi 


i 


- *- - 




(12 




\ 





(^2222) (Vi 311 ) 

(^331 1) 

FIG. 12. One-loop diagrams for (S 4 ) c corresponding to D2222 and D3311 (see Eqs. (4.26) and (4.28)). 



qi 



q4 

qi ' 



q2T + qsr + q2T 

'< / q3 \ ! q4 ! q3 

► • • ► ► • ► • ► ► • 



(^3221) (^3221) (^3221) 

FIG. 13. One-loop diagrams for (S 4 ) c corresponding to D3221 (see Eq. (4.27)). 



qi qi 

qi ' *~ •* 

v _ q2 _ ' . ' q2 



. \ q 3 . \q2 

q4T V + q 4 T n \ + qsr 

1 ~> 1 q3 n 1 



q4 



(fsiii) (V{ 211 ) 

(^4211) 

FIG. 14. One-loop diagrams for (S 4 ) c corresponding to D5111 and D4211 (see Eqs. (4.30) and (4.29)). 
We now turn to the calculation of cumulants for the different perturbation theory kernels derived in Section II. 
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V. APPLICATIONS: LOOP CORRECTIONS FOR DENSITY FIELD CUMULANTS 



Using the results derived in Section IV, we can now calculate loop corrections of 1-point cumulants for the exact 
dynamics and for the different non-linear approximations. All we need to specify is the linear power spectrum P\{k, r) 
corresponding to the linear density fluctuations. We take this to be a truncated power-law: 

n ^' Tj = \0 otherwise, (bA > 

where A is a constant amplitude, n is the spectral index, e is an infinitesimal infrared cutoff and k c is an ultraviolet 
cutoff. The wave- number cutoffs are introduced in order to regularize the radial integrals. The ultraviolet cutoff k c 
is required because at high k one reaches the strongly non-linear regime and perturbation theory breaks down. The 
spectral index is taken to be in the range —2 < n < 2. The lower limit n = —2 is chosen because for n = —3 there is 
an infrared (logarithmic) divergence in the expansion parameter o^. The upper limit n = 2 is taken to cover a range 
in n of physical interest. We have also considered other schemes for regularizing the integrals, such as dimensional 
regularization. 

For the exact dynamics and the Zel'dovich approximation, it turns out that the results for the scaled one-point 
cumulants S p can be written in terms of aj , independent of the cutoffs aside from the implicit cutoff-dependence of 
a 1 itself. In other words, the S p parameters for these cases are infrared-finite, and their ultraviolet divergences can 
be absorbed into the expansion parameter aj ~ k" +3 . 

A. Loop Corrections for Non-Linear Approximations 

None of the kernels in the different non-linear approximations derived in Section II C involve scalar products of 
wave- vectors in their denominators. Thus, the relevant angular integrals are: 

J dfii J di2 2 J di2 3 ( qi • q 2 ) a (q 2 • q 3 ) 6 (q 3 • qi) c , (5.2) 

for 0(<j 6 ,) corrections and: 

dfii J dn 2 J dtt 3 J dQ A (qi-q 2 ) a (q2-q3) 6 (q3-qi) c (qi-q4) d (q2-q4) e (q3-q 4 ) / , (5.3) 

for 0(af) corrections. Even though these integrals are straightforward to evaluate analytically, the number of terms to 
be integrated makes calculation by hand not feasible for p > 3 in the 1-loop case. For example, in the ZA, calculation 
of 1-loop corrections to S3 and S4 involve 150 and 2431 terms to be integrated, respectively. The situation worsens 
considerably for higher p and/or higher number of loops. Therefore, it is essential to perform the calculation with 
the aid of a Computer Algebra System. We have used Mathematica [64] to generate the different perturbation theory 
kernels according to their recursion relations, and then developed a symbolic integration routine that uses the results 
derived in Appendix A (see Eq. (A2) and Eq. (A3)). The resulting radial integrals are straightforward, since they 
only involve powers of radial components of wave- vectors and therefore decouple into a product of one-dimensional 
integrations. The final result is then expanded in a Laurent series of e/k c and the limit e — > is taken afterwards, 
except for the divergent terms. 

1. Loop Corrections for a 2 
We now proceed to calculate one-loop corrections to the variance of the density field. Let: 

Using the expressions for the perturbation theory kernels and doing the angular integrals using Eq. (A2) and Eq. (A3) 
we get (see Eqns. (4.4)): 
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and: 



15 6 qt; 6 q{ 

iP _ 1219 49_rf 49_fl| (56) 
22 " 1500 + 600 g 2 + 600 g 2 ' 

i22 " 12 + 24 g 2 + 24 g 2 ' 



o gj 



= -i - H (5.10) 



Now, integrating over radial components of wave vectors and dividing by af we get: 



1 Q 1 

V z=— + -k 2 k- 2 , (5.11) 
22 15 3 K ! 

1219 49 ^r- 



V^p = + k 2 k~ 2 , (5.12) 

22 1500 300 v ; 

7 1 T^T- 



and: 



V F 9 F = — + —k 2 k~ 2 , (5.13) 
12 1 2 v ; 



VI, = -h 2 k~ 2 , (5.14) 
1 19 



V^ p = k 2 k~ 2 , (5.15) 

31 15 105 K ! 

1 1- 



where we have defined: 



V FF = k 2 k~ 2 , (5.16) 

si 9 9 v ; 



k m = \ / d<7 2 (q) q m . (5.17) 



Summing the two diagrams we obtain the final result: 



4" = if. (5.18) 

17 1 



4 } F = ^-^ 2 k- 2 . (5.20) 



36 36 
Since: 

f (A+l)/3 + 0(A" 1 ) if n = -2, 

k 2 k~ 2 = I In A + C(A" 2 In A) ifn = -l, (5.21) 

( (n + 3) 2 /(n + l)(n + 5) + C^A""- 1 ) if n > 0, 

where A = fc c /e, from Eq. (5.11), Eq. (5.14), and Eq. (5.21), we see that for ZA the infrared divergences present in T>si 
and V22 for n = —2, —1 as e — > exactly cancel each other, leaving a finite correction which is independent of spectral 
index. However, for the LPA and FFA, infrared divergences do not cancel and the 1-loop correction s^ 1 ) diverges for 
n = — 2, — 1. We postpone a discussion of the physical reason behind cancellation of infrared divergences to Sec. VII. 
Note on the other hand that for n = 0,1,2 the corrections are finite and almost independent of spectral index. In 
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these cases, the loop corrections to a 2 are most important for the ZA, followed by the LPA, with FFA having the 
smallest corrections. This trend in fact agrees with the behavior of tree-level calculations [52,53,51], although the 
inclusion of smoothing of the fields changes this situation somewhat [51]. Tables I, II, and III summarize the results 
for 1-loop corrections to a 2 in the different non-linear approximations as a function of spectral index. 



TABLE I. One-loop corrections to a 2 in the Zel'dovich Approximation (A = k c /e) 



n 


Z>31 


V22 


= V31 + V22 


-2 


-(1/9)A - 1/9 


(1/9)A + 62/45 


19/15 « 1.266 


-1 


-(1/3) In A 


(1/3) In A + 19/15 


19/15 « 1.266 





-3/5 


28/15 


19/15 Rj 1.266 


1 


-4/9 


77/45 


19/15 rj 1.266 


2 


-25/63 


524/315 


19/15 k, 1.266 



TABLE II. One-loop corrections to a 2 in the Linear Potential Approximation (A = k c /e) 



_n V31 V22 = V31 +V22 

-2 -(19/315)A - 8/63 (49/900)A + 1951/2250 -(37/6300)A + 11657/15750 

-1 -(19/105) In A- 1/15 (49/300) In A + 1219/1500 -(37/2100) In A + 373/500 

-206/525 83/75 5/7 r; 0.714 

1 -97/315 4637/4500 22759/31500 rj 0.722 

2 -622/2205 1133/1125 39967/55125 rj 0.725 



TABLE III. One-loop corrections to a 2 in the Frozen Flow Approximation (A = k c /e) 



_n V31 V22 = V31 +V22 

-2 -(1/27)A - 4/27 (1/36)A + 11/18 -(1/108)A + 25/54 

-1 -(1/9) In A - 1/9 (1/12) In A + 7/12 -(1/36) In A + 17/36 

-14/45 11/15 19/45 rj 0.422 

1 -7/27 25/36 47/108 rj 0.435 

2 -46/189 43/63 83/189 rj 0.439 
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For ZA, we also compute two and three loop corrections to the variance. The calculations follow along the same 
lines as shown above for the one-loop case. Divergences from different diagrams present for n = — 2, — 1 also cancel in 
this case, leaving a finite correction independent of spectral index. These results are shown in tables IV and V as a 
function of spectral index. 



TABLE IV. Two-loop corrections to a 2 in the Zel'dovich Approximation (A = k c /e) 



n 


£>24 


u 3,3 


^33 - ^51 




2) 


-2 
-1 

1 

2 


-(1/45)A 2 - (26/45)A - 26/45 
— (4/27) In 2 A — (5/3) In A 
-24/7 
-22/9 
-11000/5103 


(1/90)A 2 + (17/30)A + 233/90 
(2/27) In 2 A + (5/3) In A + 91/45 
3299/630 
196/45 
208819/51030 


(1/180)A 2 + (1/180)A + 1/180 
(1/27) In 2 A 
3/28 
1/18 
125/2916 


91/45 r 
91/45 r 
91/45 r 
91/45 r 
91/45 r 


d 2.022 
d 2.022 
d 2.022 
d 2.022 
d 2.022 



TABLE V. Three-loop corrections to a 2 in the Zel'dovich Approximation (A = k c /e) 



n 


2?62 — 2?44 


Z>71 


-2 
-1 

1 

2 


(1/756)A 3 + (27/350)A 2 + (289/2700)A + 647/4725 
(1/54) In 3 A + (6523/11340) In 2 A 
4832/2625 
914/945 
6231200/8251551 


-(A 3 + A 2 + A + l)/4536 
-(1/324) In 3 A 
-1/72 
-2/405 
-625/192456 


n 




u b3 


-2 
-1 

1 

2 


(1/1134)A 3 + (2162/14175)A 2 + (2791/1134)A + 65404/10125 
(1/81) In 3 A + (6523/5670) In 2 A + (506/75) In A + 1477/375 
154768/7875 
1051013/70875 
4625937328/343814625 


-(A 3 + A 2 + A + 1)/1512 
-(1/108) In 3 A 
-1/24 
-2/135 
-625/64152 


n 


u b3 


S < 3 > 


-2 
-1 

1 

2 


-(1/378)A 3 - (413/1350)A 2 - (25273/9450)A - 8801/3150 
-(1/27) In 3 A - (6523/2835) In 2 A - (506/75) In A 
-101537/5250 
-20168/1575 
-181755199/16503102 


1477/375 Rd 3.939 
1477/375 Rd 3.939 
1477/375 Rd 3.939 
1477/375 Rd 3.939 
1477/375 Rd 3.939 



al + 3.939 <rf + O{<Tl ), (5.22) 

(5.23) 

1), and 

(5.24) 

1). These results are displayed in Fig. 



We can summarize the results for the variance by writing: 

a 2 = aj + f h at + + ^af + 0{<r\°) « aj + 1.266 + 2.022 

for the ZA, 

a 2 ^aj + 0.720 <r £ 4 + 0(<r £ 6 ), 

for LPA (within 1% approximately independent of spectral index when n > — 

a 2 Rd aj + 0.431 + 0(a 6 i ), 

for FFA (within 2% approximately independent of spectral index when n > — 
15. 
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FIG. 15. The ratio of the non-linear to the linear variance a 2 j ' <j\ as a function of the linear variance <j\ for the exact 
dynamics (ED) to one loop, the Zel'dovich approximation at one (ZA), two (ZA2), and three loops (ZA3), and the LPA and 
FFA approximations to one loop. 



2. Loop Corrections for S3 

Table VI presents the results of tree-level S3 for different dynamics. One-loop corrections to S3 for the non-linear 
approximations are calculated similarly to the variance case, the only complication being the increase in the number 
of terms to be integrated. The results of these calculations are presented in tables VIII, VII and IX for ZA, LPA, 
and FFA respectively. The general trend seen for corrections to a 2 is also valid here, namely, ZA has the largest 
corrections followed by LPA and then FFA. In this case, however, infrared divergences in the FFA < S 3 > c cancel 
with the divergences in the denominator (<r 4 ), making the coefficient S^ finite (see table IX). 

For the ZA we also compute two-loop corrections, shown in table X. We can summarize the results for ZA as 
follows: 

q 52 5224 

S3 = 4 + — o-J + —c-t + 0{4) « 4 + 4.693 a 2 + 13.931 af + 0{4), (5.25) 
or by inverting Eq. (5.22): 

5 3 (<t 2 ) = 4 + °^a 2 + + 0(a e ) « 4 + 4.693 a 2 + 7.986 a 4 + 0(a e ), (5.26) 

as a function of the non-linear variance. For the ZA, the 0(cr 2 ) correction S^ was previously estimated numerically 
in [20] to be S^ « 4.6, in very good agreement with Eq. (5.25). For LPA we have: 

S 3 « 3.400 + 2.124 cr'i + O(o- 4 ), (5.27) 

(n > —1), and for FFA we obtain: 

S3 = o + o- 2 + 0(at). (5.28) 
Fig. 16 gives a visual summary of these results. 
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2 U i i i I i i i I i i i I i i i I i i i L 

0.2 0.4 0.6 0.8 1 

FIG. 16. The skewness S3 as a function of the linear variance <j\ for the exact dynamics (ED) to one loop, the Zel'dovich 
approximation at one (ZA) and two loops (ZA2), and the LPA and FFA approximations to one loop. 



TABLE VI. Tree-level S3 and S4 for different dynamics in Perturbation Theory 



Dynamics S\ = V 2 n P2211 P3111 S\ = Vtua + Vzm 

ED 34/7 « 4.857 4624/147 2728/189 60712/1323f» 45.890 

ZA 4 64/3 80/9 272/9f» 30.222 

LPA 17/5 = 3.4 1156/75 1828/315 33416/1575 « 21.216 

FFA 3 12 4 16 



TABLE VII. One-loop corrections to S3 in the Linear Potential Approximation (A = k c /e) 



_H ^222 Vln 

-2 -(343/4500)A + 73351/56250 (118/225)A + 53834/7875 -(323/1575)A - 136/315 

-1 -(343/1500) In A + 17523/12500 (118/75) In A + 16568/2625 -(323/525) In A - 17/75 

3027/3125 24002/2625 -3502/2625 

1 120997/112500 66224/7875 -1649/1575 

_2 31163/28125 64454/7875 -10574/11025 

V^l -2s' 1 ' Sf 

-2 -(5353/18900)A - 5207/9450 (629/15750)A - 198169/39375 -(1/6750)A + 1254763/590625 

-1 -(5353/6300) In A- 241/900 (629/5250) In A - 6341/1250 -(1/2250) In A + 836567/393750 

-14153/7875 -34/7 418126/196875 « 2.123 

1 -26473/18900 -386903/78750 2509001/1181250 ^2.124 

2 -42313/33075 -1358878/275625 1254538/590625 « 2.124 
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TABLE VIII. One-loop corrections to S3 in the Zel'dovich Approximation (A = k c /e) 



_n V222 VI21 V\L 2>4ii Sj lh 

-2 -(2/9)A + 436/225 (4/3) A + 14 -(4/9)A - 4/9 -(2/3)A - 2/3 352/75 ^4.693 

-1 -(2/3) In A + 54/25 41nA + 38/3 -(4/3) In A -2 In A 352/75 ^4.693 

24/25 298/15 -12/5 -18/5 352/75f»4.693 

1 286/225 18 -16/9 -8/3 352/75f»4.693 

2 2152/1575 122/7 -100/63 -50/21 352/75f»4.693 

a See Eq. (4.23b). In this case 2s' 1 'S , 3°' = 152/15 independent of spectral index. 





TABLE IX. One-loop corrections to S3 


in the Frozen Flow Approximation 


(A = k c /e) 


n 


2-*222 


^321 


^321 


-2 


-(1/36)A+ 17/18 


(2/9)A + 34/9 


-(1/9)A - 4/9 


-1 


-(1/12) lnA + 35/36 


(2/3) In A + 32/9 


-(1/3) In A - 1/3 





37/45 


214/45 


-14/15 


1 


31/36 


40/9 


-7/9 


2 


55/63 


274/63 


-46/63 


n 


Z>411 




c(l) 

^3 


-2 


-(5/36)A - 1/2 


(1/18)A - 25/9 


1 


-1 


-(5/12) In A - 13/36 


(1/6) In A - 17/6 


1 





-10/9 


-38/15 


1 


1 


-11/12 


-47/18 


1 


2 


-6/7 


-166/63 


1 



TABLE X. Two-loop corrections to S3 in the Zel'dovich Approximation (A = k c /e) 



n 


^332 


^332 


^332 


-2 
-1 

1 

2 


-(8/81)A 2 - (2242/2025)A + 36898/3375 
-(106/135) In 2 A - (644/225) In A + 13606/1125 
35446/7875 
70514/10125 
34162714/4465125 


-(34/405)A 2 - (344/405)A - 41/45 
-(19/27) In 2 A - (19/9) In A 
-211/35 
-326/81 
-124025/35721 


(1/81)A 2 + (2/81)A + 1/27 
(1/9) In 2 A 
9/25 
16/81 
625/3969 


n 


u 422 


T)JJ 

^422 


Z>611 


-2 
-1 

1 

2 


(56/405)A 2 + (7016/2025)A + 52696/3375 
(154/135) In 2 A + (2182/225) In A + 13562/1125 
260696/7875 
272678/10125 
112252328/4465125 


(4/45)A 2 - (116/225)A - 52/75 
(94/135) In 2 A - (196/175) In A 
-2232/875 
-524/225 
-78632/35721 


(32/405)A 2 + (46/405)A + 4/27 
(82/135) In 2 A 
1632/875 
404/405 
4000/5103 


n 


^431 


^431 


V 111 

^431 


-2 
-1 

1 

2 


(8/45)A 2 + (1754/225)A + 2408/75 
(58/45) In 2 A + (568/25) In A + 364/15 
181204/2625 
12736/225 
3149996/59535 


-(10/81)A 2 - (1219/405)A - 136/45 (14/135)A 2 + (19/135)A + 8/45 
— (23/27) In 2 A — (388/45) In A (7/9) In 2 A 
-3154/175 414/175 
-5186/405 34/27 
-402910/35721 11750/11907 


n 


V 1 

u b2\ 


^521 


c(2)a 
°3 


-2 
-1 

1 

2 


-(128/405)A 2 - (11861/2025)A - 1349/225 
-(328/135) In 2 A - (3647/225) In A 
-32057/875 
-51844/2025 
-114469/5103 


(1/45)A 2 + (1/45)A + 1/45 
(4/27) In 2 A 
3/7 
2/9 
125/729 


5224/375 « 13.931 
5224/375 « 13.931 
5224/375 « 13.931 
5224/375 « 13.931 
5224/375 « 13.931 


a See 


Eq. (4.24). In this case 2s (1) sf + ((s (1) ) 2 +2s (2) ) 


Sf = 12932/375 independent of spec 


:tral index. 
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3. Loop Corrections for S4 



Table VI summarizes the tree-level S 4 values for the different dynamics considered here. Results of one-loop correc- 
tions to S 4 for the non-linear approximations are shown in tables XI, XII, and XIII for LPA, ZA, and FFA respectively. 
Again, the same comments apply in this case about the size of loop corrections in the different approximations and 
the cancellation of infrared divergences. We can summarize the results for ZA as follows: 

272 1 1 0822 

S 4 = — + o-\ + O{o-f) Rj 30.222 + 98.508 aj + 0{af). (5.29) 

9 1 125 

For LPA and FFA we have, respectively ( n > —1): 

S 4 « 21.216 + 37.120 aj + 0{af), (5.30) 
S 4 « 16+ 15.00 aj + O(o-f). (5.31) 

For the ZA, was previously estimated numerically in [20] to be « 100, in very good agreement with Eq. (5.29). 
The S 4 results are displayed graphically in Fig. 17. 



TABLE XL One-loop corrections to S4 in the Linear Potential Approximation (A = k c /e) 





7X1 1 1 

^5111 


^4211 


-2 
-1 

1 

2 


-(761962/779625)A - 197752/111375 
-(761962/259875) In A - 207434/259875 
-7894828/1299375 

-13346/2835 
— 234051 64/5457375 


(4163/2250)A + 2384653/118125 
(4163/750) In A + 4332191/236250 
3346306/118125 
6080651/236250 
2946658/1 1 81 25 

Li U t: \J\J tj O / _L_LO_LZ;tJ 




T) 1 

■^3221 


■nil 
^4211 


-2 
-1 

n 
u 


-(133/10125)A 2 - (13937/10125)A + 11666324/590625 
-(133/1125) In 2 A - (1519/375) In A + 4162358/196875 

9^1 zlQ9 / I 0^87^ 


-(182002/70875)A - 354076/70785 
-(182002/23625) In A - 8194/3375 

1 Q9zl&n& / I 1 81 9^ 

— lyziouo/ iioizo 


1 
1 

2 


97^1Q199/l 771 87^ 

98fi1 QQ7? /I 771 8.7^ 


— yuuuoz / / Uo / o 
— ^7^4^(SS/4Q(S1 9^ 


n 


■^3221 


-pin 

■^3221 


-2 
-1 

1 


(8024/3375)A + 3660712/118125 
(8024/1125)A + 1126624/39375 
1632136/39375 
4503232/118125 


-(21964/23625)A - 9248/4725 
-(21964/7875) In A - 1156/1125 
-238136/39375 
-112132/23625 


2 


4382872/118125 


-719032/165375 


n 


■^3311 


■^3311 


-2 
-1 

1 

2 


(722/99225)A 2 + (6568/3969)A + 425476/23625 
(722/11025) In 2 A + (6028/1225) In A + 901154/55125 
7005592/275625 
11423266/496125 
542076664/24310125 


-(17366/33075)A - 7312/6615 
-(17366/11025) In A - 914/1575 
-188284/55125 
-88658/33075 
-568508/231525 


n 


2-*2222 




-2 
-1 

1 

2 


(2401/405000) A 2 + (2401/20250)A + 43262149/8437500 
(2401/45000) In 2 A + (2401/7500) In A + 28141141/5625000 
4044362/703125 
279681269/50625000 
34549018/6328125 


(309098/275625)A - 194765156/4134375 
(309098/275625) In A - 3116042/65625 
-33416/735 
-190128686/4134375 
-1335537272/28940625 


n 


c(l) 




-2 
-1 



(1369/19845000)A 2 - (43621/27286875)A + 56278351337/1515937500 
(1369/2205000) In 2 A - (63181/12127500) In A + 112561340299/3031875000 
14067375968/378984375 « 37.119 




1 


92081148881/2480625000 « 37.120 




2 


6204054381248/167132109375 « 37.121 
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TABLE XII. One-loop corrections to S4 in the Zel'dovich Approximation (A = k c /e) 

" X>5111 ^4211 ^4211 ^3221 

~2 -(28/9)A - 28/9 (56/9)A + 512/9 -(64/9)A - 64/9 -(4/81)A 2 - (404/81)A - 5152/135 

-1 -(28/3) In A (56/3) In A + 152/3 -(64/3)lnA -(4/9) In 2 A - (44/3) In A + 216/5 

-84/5 1264/15 -192/5 384/25 

1 -112/9 680/9 -256/9 9256/405 

2 -100/9 656/9 -1600/63 498304/19845 

V n V 111 V* V 11 

" ^3221 ^3221 ^3311 ^3311 

~l (64/9)A + 224/3 -(64/27)A - 64/27 (2/81)A 2 + (424/81)A + 1287/27 -(40/27)A - 40/27 

-1 (64/3) In A + 608/9 -(64/9)A (2/9) In 2 A + (140/9) In A + 380/9 -(40/9) In A 

4768/45 -64/5 15962/225 -8 

1 96 -256/27 5132/81 -160/27 

2 1952/21 -1600/189 242330/3969 -1000/189 

n V^T 2 S [ l u 

^2 (2/81)A 2 + (40/81)A + 34516/3375 110822/1125 « 98.508 

-1 (2/9) In 2 A + (4/3) In A + 10922/1125 110822/1125f» 98.508 

14432/1125 110822/1125 « 98.508 

1 120298/10125 110822/1125f» 98.508 

2 5760352/496125 110822/1125f» 98.508 



a See Eq. (4.32b). In this case 3s' 1 ' 5^°' = 5168/45 independent of spectral index. 

TABLE XIII. One-loop corrections to 54 in the Frozen-Flow Approximation (A = k c /e) 

" X>5111 ^4211 ^4211 ^3221 

~ -(142/405)A - 92/81 (97/162)A + 635/81 -(10/9)A - 4 -(1/243)A 2 - (101/243)A + 908/81 

-1 -(142/135) In A- 106/135 (97/54) In A + 391/54 -(10/3) In A - 26/9 -(1/27) In 2 A - (11/9) In A + 314/27 

-1808/675 1414/135 -80/9 6284/675 

1 -886/405 1561/162 -22/3 2414/243 

2 -5776/2835 5318/567 -48/7 120524/11907 

n V 11 V 111 V 11 IT 1 

" ^3221 ^3221 ^3311 ^3311 

~ (8/9)A + 136/9 -(4/9)A - 16/9 -(2/9)A - 8/9 (2/729)A 2 + (424/729)A + 1900/243 

-1 (8/3) In A + 128/9 -(4/3)A - 4/3 -(2/3) In A - 2/3 (2/81) In 2 A + (140/81) In A + 586/81 

856/45 -56/15 -28/15 21112/2025 

1 160/9 -28/9 -14/9 6986/729 

2 1096/63 -184/63 -92/63 333176/35721 

~ (1/648)A 2 + (5/162)A + 329/108 (4/9)A - 200/9 (1/5832)A 2 + (2/3645)A + 14581/972 

-1 (1/72) In 2 A +(1/12) In A + 217/72 (4/3) In A - 68/3 (1/648) In 2 A + (1/1620) In A + 48601/3240 

722/225 -304/15 30388/2025 « 15.006 

1 2041/648 -188/9 437513/29160 « 15.004 

2 12434/3969 -1328/63 2679652/178605 « 15.003 
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FIG. 17. The kurtosis 54 as a function of the linear variance <j\ for the exact dynamics (ED) at tree level, the Zel'dovich 
approximation at one loop (ZA), and the LPA and FFA approximations to one loop. 



B. Loop Corrections for the Exact Dynamics 



The calculation of loop corrections in this case is technically much more complicated than in the non-linear ap- 
proximations, due to the appearance of logarithmic terms coming from solving the exact Poisson equation. Since 
the Green's function of the Poisson equation is non-local, this introduces a weak dependence of the loop corrections 
on the spectral index n. This is the fundamental difference between exact gravitational instability and non-linear 
approximations [65,66]. The principles of calculation, however, are exactly the same, and we refer the reader to 
Appendices A and B for the details. Infrared divergences present in individual diagrams for n = —2, —1 are canceled 
when the sum over diagrams contributing to a given correction is done. The mechanism behind this cancellation is 
discussed in Sec. VII. 



1. Loop Corrections for a 2 

The results of angular integration of Eqs. (4.4a) and (4.4b) in this case are: 

_ 1219 \q\ \q\ 
/22 "73^ + 6^f + 6^' 



(5.32) 



j = 25 1 g| lq 2 79 rj (g? - g 2 2 ) 3 (7g? + 2g 2 2 ) | gl + g2 | 
31 126 42 gf 12 q\ 252 q 2 168^^1 ° |gi - g 2 |' 

Now, integrating over radial components of wave vectors and dividing by af we get: 



where: 



— 1219 l-pj-r^j 
735 3 
25 25^r- 



1 v-r- 



P31 = k 2 k~ 2 + —k 4 k~ 4 ln(Ar, fc" 1 ) + — ln(P, k~ 3 ) ln(k 5 , k~ 5 ) 

126 63 42 42 v ' ; 28 v ; 84 v ' ; 



ln(k°,k<>) = \( do-f^l) [ do- 2 M 2 ) g? q\ ln^i 
°~l J J |gl 



gl + g2| 



g2| 



(5.33) 

(5.34) 
(5.35) 

(5.36) 
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In this case the cancellation of infrared divergences can be readily checked by expanding Eq. (5.33) for small q\. 



116 \q\ 188 q\ 
315 ~ 3 qf ~ 735 ql 



2 7,^2 ' 0{ qi /q 2 )\ (5.37) 



with I31 being convergent as q 2 — ?> (in fact, I31 oc q\ as q 2 — > 0, see Eq. (4.4b) and property (1) in Sec. II B). In 
Table XIV we show the results for individual diagrams and their sum. We see that the logarithmic terms in I31 are 
responsible for the very weak n-dependence in the final results. This effect does not happen at the tree-level, because 
tree diagrams correspond to taking the angular average of the perturbation theory kernels, which results in constants 
independent of the g 8 -'s [19] and therefore independent of the amount of power at a given scale. Recall that the loop 
corrections in the ZA were independent of n. This was due to the absence of logarithmic terms, which come from 
the angular dependence in the denominator of perturbation theory kernels and appear only when solving the exact 

Poisson equation. 

From Eqns. (5.17) and (5.36) it follows that k m k~ m and \n(k m ,k~ m ) are invariant under n — » —(n + 6). This 
implies that for the exact dynamics the results for the 1-loop variance are symmetric about n = — 3. The large n 
asymptotes can be determined from the fact that: 

[n + 3) 2 — m l 



for n + ±m. Therefore, k m k m becomes unity as ±n increases. Similarly, using the expansion: 



In 



kl + 92| _ 2 f q < 



ki - 92 1 ^ 2fc + 1 \ q> 



2k + l 



where q < (g>) is the smaller (larger) of { q\, q 2 }, we obtain: 

00 „ 

ln{k m ,k~ m ) = ^ j- k m +' 2k + 1 k- m ~ 2k -\ (5.40) 

k = 

and therefore the logarithmic terms in Eq. (5.33) cancel each other in the large ±n limit. This leaves us with the one 
loop result: 

4007 

v 1 = °l + ^ 4 + 0(4) « aj + 1.817 at + O(af), (5.41) 

in the large ±n limit. 

We can therefore summarize the results for the one-loop corrections to the variance: 

a 2 = <t\ + 1.82 at + 0(a 6 i ), (5.42) 

which is independent of n to better than 1%. This analytic result was recently confirmed by numerical integrations 
given in [32]. Comparing the ED result to the equivalent corrections obtained in the non-linear approximations, we 
see that the latter underestimate the exact results by 30 % (ZA), 60 % (LPA), and 76 % (FFA) (see Fig. 15). 



TABLE XIV. One-loop corrections to a 2 in the Exact Dynamics (A = k c /e) 



n 


Z>31 


V22 


= V31 +V22 


-2 
-1 

1 

2 


-(1/9)A + (17/189 - tt 2 /336) 
-(1/3) In A + (1634/6615 - (In 2)256/2205) 
-179/315 +3tt 2 /224 
-7726/19845 + (In 2)1024/6615 
-17/105 - 5tt 2 /672 


(1/9)A + 3902/2205 
(1/3) In A + 1219/735 
332/147 
4637/2205 
4532/2205 


12301/6615 - tt 2 /336 ~ 1.830 
2521/1323 - (In 2)256/2205 « 1.825 
3727/2205 + 3?r 2 /224 « 1.822 
34007/19845 + (In 2)1024/6615 « 1.821 
835/441 - 5tt 2 /672 « 1.820 
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2. Loop Corrections for S3 



The one-loop corrections to S3 are calculated using the expressions obtained in Sec. IV C and the results of integrals 
given in the Appendices. Although the following results are much more complex than in the case of the variance, 
their structure is similar, in the sense that logarithmic terms introduce a small variation of with n and infrared 
divergences in individual diagrams (for n = — 2, — 1) also cancel when summing all the contributions. For n = — 2 we 
have: 



222 



-A- 



293404 
77175 ' 



(5.43) 



^7— 32 , 



3^ 
112 



ln(A) 



25223 
1080 



1470 16 



C(3), 



(5.44) 



34 



VlL = -— A- 

321 g 3 



578 



17^ 



1323 1176 



(5.45) 



^ — 16 , 3tt 2 , , 4 . 
= ~ 21 ~ 112 ^ 



4951280852713 1035171259tt 2 



1 



■ln 3 (2) 



6571901952000 

50392659319 , , n . 

ln(3) 



25035816960 
7577tt 2 



1617 v ' 10709766144000 
7577 , 3 275115257 

"37255680 (S} 
21 9 

:Li 3 (l/2)H 



ln(3) 



34496 



10431590400 
2231753 CIS), 



37255680 
1 3 



9tt 4 



1576973 



ln(2) 



32^ 2 



68992 20697600 

275115257 , , n , , , , 

+ 10431590400 ln(2)ln(3) 
7577 

^) + 6209280 Ll3 (- 1 / 3 )- 



ln(2) 



544 



8085 v ' 40425 
275115257 2 

20863180800 n 1 ' 
7577 . 

3104640 Ll3(1/3) 



In 2 (2) 



12418560 



(5.46) 



(0) 



836468 17tt 2 



(5.47) 



46305 588 ' 

where Li n (x) denotes the polylogarithm of order n (see Eq. (B8) ) and C( x ) is l ne Riemann zeta function. Therefore, 



S 3 1} (-2) 



67574960100073 690290107tt 2 



9tt 4 



6571901952000 
50392659319 

+ 10709766144000 
7577 



25035816960 68992 
7577tt 2 



ln(3) 



ln(3) 



1576973 
20697600 
275115257 



ln(2) 



32^ 2 
8085 



ln(2) 



544 m 2 ( 2 ) 



1 



ln 3 (3) 



37255680 
10.034 



37255680 v ' 10431590400 
275115257 13 

Ll 2(l/3) " 77^ Ll 3(-l/ 2 ) 



ln(2)ln(3) 



10431590400 
^6727 

12418560 sv ; 



4928 



7577 
6209280 



40425 
275115257 
20863180800 

Li 3 (-l/3) - 



1617 



ln 3 (2) 



In 2 (3) 

7577 
3104640 



Li 3 (l/3) 



(5.48) 



For n 



-1 we have: 



V, 



ln(A) 



34506 
8575 ' 



(5.49) 



32 



2>£ 21 = y ln(A) 



5010232 68512 



231525 77175 



ln(2), 



^321 — 



55556 8704 
' 46305 ~ 15435 ^ >' 



(5.50) 
(5.51) 



- — 16, , 1N 286837213 
1,411 = -- ln(A)+ 142619400 



324178 
1037575 



73963tt 2 



1120128 



95079600 660275 



■■•< 2 >-li> 2 >+S» 2 < 2 >-^ 2 > 



ln(3) 



8 



8085 



ln(2)ln(3) 



8085 



In 2 (3) 



8 18 249 

— Li 2 (l/3) + -Li3(l/2)-— C( 3)) (5.52) 



31 



(1) (o) __ 171428 17408 
18 ^ ~ 9261 + 15435 (j ' ( ^ 3j 



and therefore: 



n(1) , 1478165197 73963tt 2 161632, , n . 27tt 2 , , o . 19968, 2 , n , 3 , 3 , n , 324178, 
Si' (-1) = ln(2) ln(2) + lri (2) ln 3 (2) ln(3) 



8 , , 4 , 2 8 18 T . _ 249 



■ ln(2) ln(3) - —— ln 2 (3) - — — Li 2 (l/3) + — Li 3 (l/2) - t^tC(3) 



8085 w w 8085 v ' 8085 w y 539 w ' 8624 
9.815 (5.54) 



For n = we have: 



24216 



1,222 = ^575-' ^ 



68608 _5W 

^ 321 " 2205 7840 ' ^ ] 

6086 51tt 2 
"2205 + 784 



2? 32i = -^ + : ^r. (5-57) 



- — _ 350978055455171 105842425tt 2 513tt 4 54676 799tt 2 _W2_ 2 131 3 

411 ~ ~ 66266678016000 + 818107136 + 1103872 ~ 471625 ^ ' + 150920 ^ ^ ~ 8575 n ^ '~ 37730 n ^ ^ 
8360129690981 , _ 893tt 2 , _ 369823547 , , _ 369823547 , 2 893 , 3 



' 35996713984000 ln(3) + 7589120 ln(3) 2191358400 ln(2) ln(3) + 4382716800 ln (3) 7589120 ln (3) 

369823547 T . , n , 2679 T . , , , n , 2679 . , n , 393 . , n , 44817 A/oN 

+ Li 2 (l/3) W-l/3) + Li 3 (l/3) + Li 3 (l/2) C(3), (5.58) 

2191358400 w ; 3794560 1 1 1 1897280 w ; 18865 w ; 7589120 s1 ; ' 1 ; 

2s (i)c(o)_ 171428 17408 

" 2S 53 --^26T + 15435 M2) ' (5 - 59) 

and therefore: 

(i) 627058864885309 20192731tt 2 513tt 4 54676 799tt 2 _W2_ 2 131 3 
3 ^ ^ ~ 66266678016000 ~ 4090535680 + 1103872 ~ 471625 ^ ' + 150920 ^ ^ ~ 8575 n ^ '~ 37730 n ^ ^ 

, 8360129690981 + _ mmw + 36982354 _ ^3 ^ 



35996713984000 v ' 7589120 v y 2191358400 w w 4382716800 v y 7589120 
369823547 T . , n , 2679 . . , , n , 2679 . , n , 393 . , n , 44817 A/oN 

+ 2191358400 Ll2(1/3) " 3794560 Ll3( " 1/3) + 1897280 Ll3(1/3) + 18865 Ll3(1/2) " 7589120 C(3) 

9.699 (5.60) 



For n = 1 we have: 



tT~ 241954 (tfiU 

1,222 = 7717^' (5 ' 61) 



^ _ 19963372 _ 153088 

^ 321_ 694575 231525 ln(2j ' (5 ' b2j 

262684 34816 
138915 + 46305 



^321 = -^^ + ^^ln(2), (5.63) 



- — 90306554812 475876271tt 2 13563074048 , 1286tt 2 , . . 1859584 , 2 . . 4 >/oN 

X» 411 = In 2 + In 2 + In 2 2 C 3 

26473726125 22187503800 344158439625 w 1867635 w 28014525 w 231 w 

52 ' m 3 (2) + ln(3) - ln(2) ln(3) + ln 2 (3) + ^Li 2 (l/3) + ^L l3 (l/2), 



266805 v ' 52026975 v ' 31185 w w 31185 v ' 31185 v ' ' 88935 

(5.64) 



32 



-2sWsi 0) 



2312476 69632 



138915 46305 



ln(2), 



(5.65) 



and therefore: 



ln(2) 



499893398528 

26473726125 ~ 22187503800 ~ 344158439625 

52 , 3 . o . 43036964, , o . 64 , , o . , , o . 
In (2) + rnnnrnrrr ln(3) - oiior ln(2) ln(3) 



1286tH 
1867635 



ln(2) 



1859584 1^(2)-A C ( 3 ) 
28014525 v ; 231 v ; 



32 



266805 
9.635 



52026975 



31185 



31185 



In 2 (3) 



64 



31185 



Li 2 (l/3) 



104 
88935 



For n = 2 we have: 



V, 



249304 



T) 1 



301492 1061tt 2 



11025 



^321 



v 411 



667518526399 909575099tt 2 276209456 



707011905600 
17807 , 2 . 



132132 



In (3) 



8701684992 
178C7 :L, 2 (l/3) 



66066 



509864355 
8 ° 55 ;C(3) 



578 

"735 ~ 
ln(2) 



37632 
85tt 2 



2352' 

19808 
~ 3776773 



In 2 (2) 



161239541 
881760880 



ln(3) 



17807 
66066 



68992 



-2sWsi 0) 



56780 85tt 2 



3087 



and therefore: 

7390246004417 21860303tt 2 



4 1) (2) 



276209456 



707011905600 



17807 
"132132 
9.561 



ln 2 (3) 



543855312 509864355 
17807 :L l2 (l/3) + ^C(3) 



ln(2) 



1176' 



19808 
3776773 



In (2) 



161239541 
881760880 



ln(3) 



17807 
66066 



Li 3 (l/2) 
(5.66) 

(5.67) 

(5.68) 

(5.69) 

ln(2)ln(3) 
(5.70) 

(5.71) 
ln(2)ln(3) 



66066 



68992 



(5.72) 

We can summarize the results for the one loop corrections to the skewness factor in the exact dynamics writing: 

S 3 « 4.86 + 9.80 crj + Oicrf), (5.73) 

which is independent of n to better than 3% when —2 < n < 2. Comparing to the equivalent corrections given in the 
non-linear approximations, we see that these underestimate the exact results by 52 % (ZA), 78 % (LPA), and 90 % 
(FFA) (see Fig. 16). 



VI. APPLICATIONS: LOOP CORRECTIONS FOR VELOCITY FIELD DIVERGENCE CUMULANTS 

In order to simplify the expressions for the moments of the divergence of the velocity field it is convenient to 
introduce the normalized velocity divergence or expansion scalar [6] given by: 

e«^, s ^i, « 6 ,) 

so that to tree-level the variance of is simply af (see Eq. (2.6)). The same expressions derived in Section IV for the 
density field are valid for the normalized velocity divergence upon replacing the kernels F n 's by G n 's. Loop expansion 
coefficients for one-point cumulants in this case follow analogous notation, namely, gives the one-loop correction 
to <t 2 , T 3 (1) the one loop correction to T3 and so on. Also, we shall refer to the amplitude of diagrams given by the 
diagrammatic rules for by Ti...j, in similar notation to T>i j in the case of the density field. We now proceed to 
present the results of loop corrections for the normalized velocity divergence. 
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A. Loop Corrections for Non-Linear Approximations 



1. Loop Corrections for a% 



One-loop corrections to the variance of the normalized velocity divergence are presented in Table XV (ZA) and 
Table XVI (LPA). Since the FFA assumes that the velocity field remains linear, all the loop corrections to one-point 
cumulants of the normalized velocity divergence vanish in this approximation. For the ZA, tables XVII and XVIII 
present results for two and three loop corrections to cr@ respectively. We can summarize these results for ZA by 
writing: 



11 



13 6 4733 8 



o-f + —o-\ + —o-l + ^r^rro-f + 0(<t} u ) Rj af + 0.733 a\ + 0.866 a\ + 1.402 af + 0(<t} u ), 



15 1 15 1 3375 
or, by inverting Eq. (5.22): 

a 2 & (a 2 ) = a 2 - ^-a 4 + ^-a* - —a* + 0(a w ) « a 2 - 0.533 ^ + 0.195 a« - 0.266 a* + 0(a 

_L ZZO _L 

as a function of the non-linear variance. Similarly, for LPA we have: 



10\ 



(6.2) 



(6.3) 



(6.4) 



which is approximately valid when n > — 1. For n = —2,-1, the one-loop correction is infrared divergent, as it was 
in the case of the variance of the density field. These results are shown in Fig. 18. 




FIG. 18. Same as Fig. 15 but for the variance of the velocity divergence. 
TABLE XV. One-loop corrections to a% in the Zel'dovich Approximation (A = k c /e) 



n 


Tsi 


722 


t {1) = Tsi +T22 


-2 


-(1/9)A - 1/9 


(1/9)A + 38/45 


11/15 « 0.733 


-1 


-(1/3) In A 


(1/3) In A + 11/15 


11/15 Rj 0.733 





-3/5 


4/3 


11/15 « 0.733 


1 


-4/9 


53/45 


11/15 Rj 0.733 


2 


-25/63 


356/315 


11/15 Rj 0.733 
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TABLE XVI. One-loop corrections to a% in the Linear Potential Approximation (A = k c /e) 



n 


Tsi 


722 


t {1) = Tsi +T22 


-2 


-(8/315)A - 8/315 


(4/225)A + 152/1125 


-(4/525)A + 96/875 


-1 


-(8/105) In A 


(4/75) In A + 44/375 


(4/175) In A + 44/375 





-24/175 


16/75 


8/105 « 0.076 


1 


-32/315 


212/1125 


76/875 « 0.087 


2 


-40/441 


1424/7875 


552/6125 Rj 0.090 



TABLE XVII. Two-loop corrections to 0% in the Zel'dovich Approximation (A = k c /e) 

n 757 7~j W=7^ t (2) 

^2 -(1/45)A 2 - (2/5)A - 2/5 (1/90)A 2 + (7/18)A + 113/90 (1/180)A 2 + (1/180)A + 1/180 13/15 rj 0.866 

-1 —(4/27) In 2 A — (17/15) In A (2/27) In 2 A + (17/15) In A + 13/15 (l/27)ln 2 A 13/15 rj 0.866 

-432/175 3277/1050 3/28 13/15 r; 0.866 

1 -26/15 112/45 1/18 13/15 r; 0.866 

2 -7760/5103 117451/51030 125/2916 13/15 r; 0.866 



TABLE XVIII. Three-loop corrections to 0% in the Zel'dovich Approximation (A = k c /e) 



n 


762 — 


T71 


-2 
-1 

1 

2 


(l/756)A a + (179/3150)A 2 + (493/6300)A + 157/1575 
(1/54) In 3 A + (4763/11340) In 2 A 
512/375 
674/945 
4585600/8251551 


-(A 3 + A 2 + A + l)/4536 
-(1/324) In 3 A 
-1/72 
-2/405 
-625/192456 


n 


'44 


'53 


-2 
-1 

1 

2 


(1/1134)A 3 + (1586/14175)A 2 + (38239/28350)A + 198028/70875 
(1/81) In 3 A + (4763/5670) In 2 A + (806/225) In A + 4733/3375 
35336/3375 
536213/70875 
6964713704/1031443875 


-(A 3 + A 2 + A + 1)/1512 
-(1/108) In 3 A 
-1/24 
-2/135 
-625/64152 


n 


'53 


t (3) 


-2 
-1 

1 

2 


-(1/378)A 3 - (2123/9450)A 2 - (677/450)A - 5009/3150 
-(1/27) In 3 A - (4763/2835) In 2 A - (806/225) In A 
-8807/750 
-35768/4725 
-106419983/16503102 


4733/3375 rj 1.402 
4733/3375 rj 1.402 
4733/3375 rj 1.402 
4733/3375 rj 1.402 
4733/3375 rj 1.402 
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2. Loop Corrections for T3 



Table XIX shows the results of tree-level calculations of X3 for different dynamics in NLCPT. One-loop corrections 
to T3 are shown in Table XX for the ZA and Table XXI for the LPA, and Table XXII presents the results of two-loop 
corrections to X3 for ZA. We can summarize the results for ZA as follows: 



-Z (7 

225 



12916 



2 - 4 <-0(<T?) 



3375 



or, by inverting Eq. (5.22): 



368 5924 

225 3375 1 ; 



-2 - 1.635 aj - 3.827 aj + 0(a^ 



-2 - 1.635 a' 2 - 1.755 <r 4 + £>(<t 6 ) 



as a function of the non-linear variance. For LPA we have: 

T 3 s» -0.9-0.19 aj + 0(af), 
approximately independent of n when n > — 1. Fig. 19 shows these results. 





-4 



-8 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 




LPA (n>-l) 


-""^^^ED ( n =2) 


ZA ^ 




N N ZA2 


ED (n=-2) ^ 




1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 



0.2 0.4 0.6 0.8 1 

FIG. 19. Same as Fig. 16 but for the skewness of the velocity divergence. 



(6.5) 



(6.6) 



(6.7) 



TABLE XIX. Tree-level T3 and T4 for different dynamics in Perturbation Theory 



Dynamics 



Tl 



(0) 



T21 



T 25 



T, 



(0) 



72211 + 731 



ED 
ZA 
LPA 
FFA 



-26/7 « -3.714 2704/147 568/63 12088/441 « 27.410 

-2 16/3 8/3 8 

-4/5 = -0.9 64/75 64/105 256/175 « 1.463 
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TABLE XX. One-loop corrections to T3 in the Zel'dovich Approximation (A = k c /e) 



n 


7222 


T 1 

'321 


T n 

'321 


Tin 


T (l)a 
1 3 


-2 


(2/9)A + 4/75 


-(8/9)A - 238/45 


(2/9)A + 2/9 


(4/9)A + 4/9 


-368/225 « -1.635 


-1 


(2/3) In A - 38/225 


-(8/3) In A - 22/5 


(2/3) In A 


(4/3) In A 


-368/225 « -1.635 





232/225 


-46/5 


6/15 


12/15 


-368/225 Rj -1.635 


1 


18/25 


-358/45 


8/9 


16/9 


-368/225 rj -1.635 


2 


328/525 


-2386/315 


50/63 


100/63 


-368/225 rj -1.635 



a See Eq. (4.23b). In this case 2t (1) T 3 (0) = -44/15 independent of spectral index. 





TABLE XXI. One-loop corrections to 


T3 in the Linear Potential Approximat 


ion (A = k c /e) 


n 


7222 


'321 


'321 


-2 


(16/1125)A + 32/9375 


-(128/1575)A - 544/1125 


(32/1575)A + 32/1575 


-1 


(16/375) In A - 304/875 


-(128/525) In A - 352/2625 


(32/525) In A 





1856/28125 


-736/875 


96/875 


1 


144/3125 


-5728/7875 


128/1575 


2 


2624/65625 


-38176/55125 


32/441 


n 


Tin 


-2t( 1 )T 3 (u) 


T (i) 

1 3 


-2 


(304/4725)A + 352/4725 


-(32/2625)A + 768/4375 


(128/23625)A - 123904/590625 


-1 


(304/1575) In A + 16/1575 


-(32/875) In A + 352/1875 


(128/7875) In A - 42368/196875 





2816/7875 


64/525 


-36608/196875 rj -0.186 


1 


1264/4725 


608/4375 


-114304/590625 rj -0.193 


2 


7936/33075 


4416/30625 


-809728/4134375 rj -0.196 



TABLE XXII. Two-loop corrections to T3 in the Zel'dovich Approximation (A = k c /e) 



n 


T 1 

'332 


T n 

'332 


■Till 
'332 


-2 
-1 

1 


(37/405)A 2 + (2626/2025)A + 1511/3375 
(97/135) In 2 A + (782/225) In A - 112/125 
6637/875 
49948/10125 


(8/135)A 2 + (46/135)A + 17/45 
(13/27) In 2 A + (11/15) In A 
99/35 
242/135 


-(1/162)A 2 - (1/81)A - 1/54 
-(1/18) In 2 A 
-9/50 
-8/81 


2 


18677123/4465125 


54185/35721 


-625/7938 


n 


r 1 

'422 


T n 

'422 


7611 


-2 
-1 


-(44/405)A 2 - (3644/2025)A - 16208/3375 
-(118/135) In 2 A - (122/25) In A - 1102/375 


-(34/405)A 2 - (292/2025)A - 124/675 
-(88/135) In 2 A - (14/225) In A 


-(8/135)A 2 - (23/270)A - 1/9 
-(41/90) In 2 A 



1 

2 


-37936/2625 
-110534/10125 
-44249264/4465125 


-1856/875 
-2348/2025 
-32896/35721 


-1224/875 
-101/135 
-1000/1701 


n 


T 1 

'431 


T n 

'431 


■Till 
'431 


-2 
-1 

1 

2 


-(52/405)A 2 - (7564/2025)A - 7228/675 
-(128/135) In 2 A - (2408/225) In A - 104/15 
-76292/2625 
-46016/2025 
-3726388/178605 


(32/405)A 2 + (529/405)A + 178/135 
(5/9) In 2 A + (164/45) In A 
1436/175 
2318/405 
59660/11907 


-(28/405)A 2 - (38/405)A - 16/135 
-(14/27) In 2 A 
-276/175 
-68/81 
-23500/35721 


n 


r 1 

'521 


T n 

'521 


T (2) a 
1 3 


-2 
-1 

1 

2 


(32/135)A 2 + (661/225)A + 2053/675 
(82/45) In 2 A + (1753/225) In A 
17167/875 
9032/675 
19777/1701 


-(1/90)A 2 - (1/90)A - 1/90 
-(2/27) In 2 A 
-3/14 
-1/9 
-125/1458 


-12916/3375 x -3.827 
-12916/3375 rj -3.827 
-12916/3375 rj -3.827 
-12916/3375 rj -3.827 
-12916/3375 rj -3.827 



See Eq. (4.24). In this case 2t (1) T 3 (1) + ((t (1) ) 2 + 2t (2) )T 3 (0) = -23426/3375 independent of spectral index. 
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3. Loop Corrections for T4 



One-loop corrections to T4 are given for the ZA in Table XXIII and for LPA in Table XXIV 
values given in Table XIX, we can summarize the results for ZA as follows: 

21 98 

T4 = 8 + — <t| + O{o-t) « 8 + 17.584 aj + 0(*t). 

For LPA we have (n > —1): 

T 4 « 1.463 + aj + 0(af). 
The kurtosis of the velocity divergence in these approximations is shown in Fig. 20. 

30 



25 
20 

H* 15 
10 

5 





0.2 0.4 0.6 0.8 1 



FIG. 20. Same as Fig. 16 but for the kurtosis of the velocity divergence. 



TABLE XXIII. One-loop corrections to T4 in the Zel'dovich Approximation (A = k c /e) 



n 


75111 


T 1 

'4211 


T n 

'4211 


T 1 

'3221 


-2 
-1 

1 

2 


-(40/27)A - 40/27 
-(40/9) In A 
-8 
-160/27 
-1000/189 


-(64/27)A - 64/27 
-(64/9) In A 
-64/5 
-256/27 
-1600/189 


(80/27)A + 1984/135 
(80/9) In A + 176/15 
416/15 
3184/135 
21088/945 


-(4/81)A 2 - (260/81)A - 832/675 
-(4/9) In 2 A - (28/3) In A + 152/75 
-1216/75 
-22696/2025 
-963904/99225 


n 


T n 

'3221 


■Till 
'3221 


T 1 

'3311 


T n 

'3311 


-2 
-1 



(64/27)A + 1904/135 
(64/9) In A + 176/15 
368/15 


-(16/27)A - 16/27 
-(16/9)A 
-16/5 


(2/81)A 2 + (184/81)A + 1498/135 
(2/9) In 2 A + (20/3) In A + 44/5 
538/25 


-(4/9)A - 4/9 
-(4/3) In A 
-12/5 


1 

2 


2864/135 
19088/945 


-64/27 
-400/189 


7324/405 
338386/19845 


-16/9 
-100/63 


n 






72222 


T (l)a 


-2 
-1 

1 






(2/81)A 2 + (40/81)A + 4756/3375 
(2/9) In 2 A + (4/3) In A + 334/375 
1504/375 
31018/10125 


2198/125 « 17.584 
2198/125 « 17.584 
2198/125 « 17.584 
2198/125 « 17.584 


2 






1385632/496125 


2198/125 « 17.584 



a See Eq. (4.32b). In this case 3t^T^ = 88/5 independent of spectral index. 



Using the tree-level 
(6.8) 



(6.9) 



"J 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 J" 




LPA (n>4) 

J 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 L 
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TABLE XXIV. One-loop corrections to T4 in the Linear Potential Approximation (A = k c /e) 



n 


75111 




T 1 

'4211 


-2 
-1 

1 

2 


-(21376/155925)A - 3712/22275 
-(21376/51975) In A - 512/17325 
-66688/86625 
-8192/14175 
-566656/1091475 




-(9728/70875)A - 11264/70875 
-(9728/23625) In A - 512/23625 
-90112/118125 
-40448/70875 
-253952/496125 


n 


T n 

'4211 




T 1 

'3221 


-2 
-1 

1 

2 


(3968/23625)A + 96512/118125 
(3968/7875) In A + 76672/118125 
183808/118125 
156032/118125 
1032704/826875 


-(1664/14175)A - 26624/590625 
-(128/7875) In 2 A - (128/375) In A + 4864/65625 
-38912/65625 
-726272/1771875 
-30844928/86821875 


n 


T n 

'3221 




■Till 
'3221 


-2 
-1 

1 

2 


(2048/23625)A + 8704/16875 
(2048/7875)A + 5632/13125 
11776/13125 
91648/118125 
610816/826875 




-(512/23625)A - 512/23625 
-(512/7875) In A 
-512/4375 
-2048/23625 
-512/6615 


n 


T 1 

' 331 1 




T n 

'3311 


-2 
-1 

1 

2 


(128/99225)A 2 + (11776/99225)A + 13696/23625 
(128/11025) In 2 A + (256/735) In A + 2816/6125 
34432/30625 
468736/496125 
21656704/24310125 




-(256/11025)A - 256/11025 
-(256/3675) In A 
-768/6125 
-1024/11025 
-256/3087 


n 


72222 






-2 
-1 

1 

2 


(32/50625)A 2 + (128/10125)A + 76096/2109375 
(32/5625) In 2 A + (64/1875) In A + 5344/234375 
24064/234375 
496288/6328125 
22170112/310078125 


(1024/30625)A - 73728/153125 
(3072/30625) In A - 11264/21875 
-2048/6125 
-58368/153125 
-423936/1071875 


n 






T (i) 


-2 
-1 

1 

2 




(32/275625)A 2 - (17408/1010625)A + 132846016/126328125 
(32/30625) In 2 A - (52928/1010625) In A + 405022048/378984375 
370578688/378984375 « 0.978 
34478368/34453125 « 1.001 
18715780352/18570234375 « 1.008 


TABLE XXV. One-loop corrections to 0% in the Exact Dynamics (A = k c /e) 


n 


T31 


T22 


t {1) = Tsi +T22 


-2 
-1 

1 

2 


-(1/9)A - (5/63 + tt 2 /112) 
-(1/3) In A + (374/2205 - (In 2)256/735) 
-113/105 +9tt 2 /224 
-5626/6615 + (In 2)1024/2205 
-83/315 - 5tt 2 /224 


(1/9)A + 2798/2205 2623/2205 - tt 2 /112 « 1.101 
(1/3) In A + 851/735 2927/2205 - (ln2)256/735 « 1.086 
1292/735 167/245 + 9?r 2 /224 « 1.078 
3533/2205 4973/6615 + (In 2)1024/2205 « 1.074 
3428/2205 949/735 - 5?r 2 /224 « 1.071 
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B. Loop Corrections for the Exact Dynamics 



We now present the results of one-loop corrections to one-point cumulants of the normalized velocity divergence for 
the exact dynamics. The same remarks made above for the density field apply here, namely, the logarithmic terms 
coming from the non-locality of the Poisson Green's function induce a small variation of the corrections with spectral 
index, and infrared divergences present for n = —2, —1 are exactly canceled when the sum over diagrams is done. 



1. Loop Corrections for a% 

Table XXV (see previous page) gives the results of one-loop corrections to the variance of for the exact dynamics. 
Comparing with the density field correction, we see that cr@ is less affected by non-linearities than a, at least in the 
weakly non-linear regime. We can summarize the results of Table XXV by: 

<7|«<7 £ 2 + l.O8<7 £ 4 + 0(<7 £ 6 ), (6.10) 

or, by inverting Eq. (5.42): 

4((T 2 ) Rj a 2 -0.74 a 4 + 0(a 6 ). (6.11) 

Comparison with the one-loop corrections in the non-linear approximations indicates that ZA underestimates the 
correct results by 32 % and LPA by 92 % (see Fig. 18). 

2. Loop Corrections for T3 

We now present the results of one-loop corrections to T3 in the exact dynamics. For n = — 2 we have: 

75Ha-§|, <•-•*) 

80 . 3tt 2 , , . , 79711 7T 2 3 



T - = "63 A " 112 " 1880" + 245 + 16« 3 >' ( ° 3) 
26 . 130 13tt 2 



40. 3?r 2 , 708960735827 810641107tt 2 9tt 4 1576973, _ 128tt 2 , 2176 



r4n 63 A+ 112 1642975488000 6258954240 + 17248 + 5174400 8085 40425 ln ^ 

A ' m 3 (2) - ln(3) + H*) + HD H*) ~ ™™1 ^(3) 



1617 v ' 2677441536000 v ' 18627840 v ' 2607897600 w w 5215795200 
73781 , 3 178450997 T . „ 13 T . , , 73781 T . , , . 73781 



111 (3) " o g ^on^nn Li 2(V3) " 77^Li 3 (-l/2) ~ o 1 n , a , n Li 3 (~ 1/3) + , K COTOn Li 3 (l/3) 



18627840 v ' 2607897600 y 1 ' 1232 ' ' ' 3104640 v 11 1552320 
219 , . . 1415461 , , , , 

+ 862i Ll3 ( 1 / 2 ) + 6209280^' ^ 

m (o) 136396 13tt 2 
- a( )T 3 0) = 15435- -W' ( ° 6) 

and therefore: 

(i) 8977705003693 628620499tt 2 9tt 4 1576973 128tt 2 2176 2 

3 ^ _ ^ ~ " 1642975488000 ~ 6258954240 + 17248 + 5174400 ^ ' ~ 8085 ^ ^ ~ 40425 n ^ ' 

4 , 3 . o . 175638489101 , , , 73781tt 2 , , o . 178450997 , , n , , , , 178450997 , 2 , , 

In (2) ln(3) H ln(3) H ln(2) ln(3) In (3) 

1617 v ; 2677441536000 v ; 18627840 v ; 2607897600 w w 5215795200 v ; 

73781 , 3/oN 178450997 T . , n , 13 T . , , , n , 73781 T . , , , n , 73781 . , n , 

In (3) Li 2 (l/3) Li 3 (-l/2) Li 3 (-l/3) + Li 3 (l/3) 

18627840 v ; 2607897600 v ' ; 1232 v ' ; 3104640 v ' ; 1552320 v ' ; 

+ ^L l3 (l/2) - ^1 C (3) 
8624 v ' ; 6209280 v ; 

« -4.643 (6.17) 
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For n = — 1 we have: 



2, 130682 



T 222 =-ln(A)-— , (6.18) 



—7- 80, , 1N 937544 19424, . . 

T - = -2T ln ( A )- 77175-- 25725 ^ 



26 , , . , 9724 6656 



To^ = — ln(A) - + ln(2), (6.20) 

321 21 1 ; 15435 5145 K h K ! 

40, , 1N 32159713 73963tt 2 5248736, 54?r 2 , . . 79872 , 12 , 3 , . 

1 = 21 ln(A) " 35654850 + 23769900 + 1980825 ln(2) + 3773 ln(2) " 660275 ln (2) + 539 ln (2) 

+ Si ln ( 3 ) -sl5M2)M3) + ^m'(3) + «-U.<l/*> - §Li 3 (l/2) + ^CO), (6-21) 

(0) = 152204 _ 13312 

3 15435 5145 w v ; 



and therefore: 



7862053 73963tt 2 1190528, , n , 54tt 2 , , o . 79872 , 2/oN 12 , 3 , o . 
Ts = " 1426194 + 23769900 + 1980825 ln(2) + 3773 M2) " 660275 ln (2) + 539 ln (2) 

+ SH " 8^5 + 8^5 ^ + ~ i^d/2) + U^ 3 ) 

« -3.577 (6.23) 



For n = we have: 



38072 



7222 = -77175' ^ 



7fr __ 14704 423^ 

/321 - 735 + 784Q > 1^5) 



2938 117tH 



—- 146302676485571 142985207tt 2 513tt 4 218704, . . 799tt 2 , . . 768 , 2 ,„. 

Taw = h h In 2 + In 2 In 2 2 

16566669504000 1022633920 275968 471625 v ; 37730 v ; 8575 v ; 

262 , 3 , o . 8360129690981, , , 893tt 2 , , o . 369823547, , n , , 369823547, 2 , , 

In (2) H ln(3) H ln(3) ln(2) ln(3) H In (3) 

18865 v ; 8999178496000 v ; 1897280 v ; 547839600 w w 1095679200 v ; 

893 , 3, n , 369823547 T . , n , 2679 . . , , n , 2679 T . , n , 1572 . , , , 

In (3) + Li 2 (l/3) W-l/3) + Li 3 (l/3) + Li 3 (l/2) 

1897280 v ; 547839600 v ' ; 948640 v ' ; 474320 w ; 18865 v ' ; 

44817 

1897280 s1 ; ' 1 ; 

m to) 8684 117tt 2 



and therefore: 



41 



T 3 (1) (0) 



43185676451389 64802527-., 
" " 16566669504000 + 1022633920 
262 3 8360129690981 
~ 18865 n ^ ' + 8999178496000 



'tt 2 513?r 4 218704, 
— + — — - - — — 1 



oidir* zio(U4 7997T 2 768 , 2 . ri . 

275968 " 471625 ln(2) + 37730 ln(2) " 8575 ln (2) 
8360129690981 893tt 2 369823547 

' 8999178496000 ^ ^ + 1897280 ^ ^ ~ 547839600 
893 3 369823547 _ 2679 

1897280 n 1 j 547839600 121 7 ' 948641 

44817 cf3) 
1897280 sv ; 
i -3.123, 



For n = 1 we have: 



168 471625 v ' 37730 v ' 8575 v ' 
893tt 2 369823547 369823547 

1897280 ^ ^ ~ 547839600 ^ ' ^ ' + 1095679200 



ln 2 (3) 



r 2 



222 



6898 
"8575' 



(6.29) 



(6.30) 



T 1 

321 



4170884 37376 



231525 77175 



ln(2), 



321 



146276 26624 



46305 15435 



ln(2), 



(6.31) 
(6.32) 



^ 150452937298 475876271tt 2 1412780667392, . . 

Till = o.^o^-or + rr ,^ r „ r „ + n , , , r n , n,wor ln ( 2 ) 



26473726125 



208 
"266805 



ln 3 (2) 



5546875950 344158439625 

172147856, 256 , , n , , , , 128 , 2 , , 

ln(3) + ln(2) ln(3) - oiior In 2 (3) 



^lln(2)- 7438336 ln 2 (2) + ii C (3) 
1867635 w 28014525 w 231 w 



52026975 



31185 



31185 



256 
31185 



Li 2 (l/3) 



416 
88935 



Li 3 (l/2), 



(6.33) 



_ 2;(1) (0) = 258596 + 53248 

3 46305 15435 v ; 



and therefore: 



116287648352 475876271tt 2 



2173099907072 

26473726125 + 5546875950 + 344158439625 
208 , 3 _ 172147856, 256 



ln(2) 



5144?r 2 
1867635 



ln(2) 



266805 
-2.871 



In (2) 



For n = 2 we have: 



52026975 



ln(3) 



31185 



1 



(6.34) 



7438336 ^ 
28014525 v ; 231 v ; 

25 Vi 2 (l/3)-^Li 3 (l/2) 



23144 
"25725' 



31185 



88935 



(6.35) 
(6.36) 



T 1 

321 



186692 247?r 2 



321 



11025 12544' 



2158 65tt 2 
2205 + 784 ' 



(6.37) 
(6.38) 



7411 



239944659679 2080959611tt 2 



176752976400 8701684992 
17807, 2 ,„, 35614, 



33033 



In (3) 



33033 



Li 2 (l/3) 



1104837824 
509864355 

17248 U h 



ln(2) 



79232 , 2 . o . 161239541 , o . 35614, , o . , , o . 
■ ln ( 2 ) + nnnAAnnnn ln(3) - oonoo ln(2) ln(3) 



3776773 



220440220 



33033 



(6.39) 



(1) (o) = 49348 _ 65^ 
3 5145 392 ' 



(6.40) 
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and therefore: 

(i) 1043821311809 74261045tt 2 1104837824 79232 2 161239541 35614 

3 ^ ~ 176752976400 + 543855312 + 509864355 ^ ' + 3776773 11 ^ ^ + 220440220 ^ ^ ~ 33033 ^ ^ 
, 17807, 2/oN 35614 T . , o . 8055 A/oN 

X ln ^) + 33033 + 33033^(V3) " T^C(3) 

« -2.579 (6.41) 

We see that the dependence of the velocity divergence skewness upon n is stronger than for the density field 
skewness. Taking n = 1 as a reference, we find that ZA underestimates the one-loop correction by 43 % and LPA by 
93 % (see Fig. 19). 



VII. CANCELLATION OF INFRARED DIVERGENCES AND GALILEAN IN VARIANCE 

As we saw in the previous sections, individual loop diagrams contain infrared divergences for spectral indices 
n = —2,-1. These divergences cancel when the total contribution to a given cumulant is computed (the diagrams 
are summed over) for the exact dynamics and the Zel'dovich approximation. However, for LPA and FFA, infrared 
divergences remain in the final answer. The issue of cancellation of infrared divergences has been considered recently 
by [67], who showed that leading infrared divergences always cancel in the contribution to the power spectrum 
to arbitrary number of loops in the exact dynamics. Although this does not prove the cancellation of subdominant 
divergences, it suggests that there is an underlying mechanism for these cancellations. In this section, we show that this 
behavior can be related to the properties of the corresponding equations of motion under a Galilean transformation. 
This is analogous to what happens in quantum field theory, where symmetries such as gauge invariance play a critical 
role in the structure and cancellation of divergences. In fact, symmetries lead to relations between correlation functions 
(the Ward Identities) which are crucial to the proof of renormalizability [68]. 

The Newtonian equations of motion of gravitational instability constitute a classical non-relativistic field theory 
and therefore are invariant under Galilean transformations (GT). In order to understand the behavior of the exact 
dynamics and the non-linear approximations under GT, we recall some basic properties of GT in comoving coordinates. 
Under a GT, the proper coordinates (r,t) transform as: 

r' = r-ut (7.1) 
V'(r') = V(r) - u, (7.2) 

where V = dr/dt, and u is a uniform velocity. In terms of comoving conformal coordinates (x, r), we have r = ax 
and V = Wx + v, where v = dx/dr, and the GT in comoving coordinates becomes: 

x' = x-uT (7.3) 

v'(x')=v(x)-u(l-F), (7.4) 

where 

T(T) = -±-r f a(T')dr'. (7.5) 
a ( T ) Jo 

Consider the transformation of individual terms in the equations of motion of the exact dynamics (see Eqs. (2.1)). 
Since for any scalar field 'f we have ^'(x') = ^(x) = \?(x' + uT) (except for the velocity field which also undergoes 
a transformation of its homogeneous mode), the partial time derivative transforms as: 

7T — ► 7T + (l-^T)u-V, (7.6) 

OT OT 

where we have used the fact that dT/dr = (1 — %T). This means that the operator 

is Galilean invariant, which in turn implies the invariance of mass conservation, Eq. (2.1a). The Galilean invariance 
of the momentum conservation equation can also be demonstrated; under a GT we have: 
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-j- — > -j- + u + 71(1 - 7LT) u, (7.8a) 
dr dr or 

Tiv — ► Tiv -71(1 -TIT) u, (7.8b) 
37-t 

V$ — ^v$-T— u, (7.8c) 

OT 

where the last expression arises because $ is the gravitational potential of the density fluctuations, related to the 
gravitational potential of the total density ip by $ = ip + (l/2)x 2 d%l dr [6]. Substituting Eq. (7.8) into the Euler 
equation, Eq. (2.1b), demonstrates the Galilean invariance of momentum conservation in a general homogeneous and 
isotropic background. For the Einstein-de Sitter model, we have T = r/3 and letting u — > 3u the GT becomes 

Wv^Wv-Wu, (7.10) 
V$ — >V$ + Wu. (7.11) 

Thus, in the Einstein-de Sitter case, the GT-dependent terms in the gravitational force and the Hubble expansion 
drag cancel each other. 

It is also useful to view the invariance properties of the equations of motion in momentum space. In Fourier space, 
the density and velocity fields transform under a GT (for the Einstein-de Sitter case) as: 

S(k) — > exp(iru • k)S(k) (7.12) 
v(k) — > exp(iru • k)v(k) - uS D (k). (7.13) 

Infrared divergences come from the v • V terms in the equations of motion [67]. These "convective" terms represent 
the time variation of fields due to the transport of fluid elements in Eulerian space. This transport is dominated by 
the homogeneous mode of the velocity field, as can be seen from the infrared behavior of the mode coupling functions 
a(k, ki) and /?(k, ki,k2) (see Eq. (2.5)). Thus these divergences are just a kinematical effect due to the fact that 
the linear fluctuations are characterized by the density field power spectrum, and therefore the rms velocity field gets 
divergent contributions from the homogeneous mode for spectral indices n = —2, —1 [67]. To see this, since we are 
interested in the infrared behavior, we can use linear theory to relate the velocity field power spectrum to the density 
power spectrum: 

<v l (k)v 1 (k')>=S D (k + k') BlJ (k) ^ -U 2 ^P l (k)8 D (k + k l ), (7.14) 

which shows that the rms velocity, v rms , becomes infrared divergent for spectral indices n = —2, —1: 

v 2 rms = J d 3 kBa(k) w -U 2 J d 3 k^^- oc J dk k n . (7.15) 

Since the infrared divergences represent a kinematic effect, they should not appear when calculating Galilean- 
invariant quantities, such as the density field (equal-time) p-point cumulants. From (7.12), under a GT these 
transform as 

<<y(ki)...<y(kp) > — ► ex P [iru-(ki + ... + kp)]<<y(ki)...<y(kp)>, (7.16) 

and are Galilean-invariant due to translation invariance, which requires X^=i k«' = 0- Similar analysis applies to 
the velocity field divergence. Non-equal-time p-point cumulants are not invariant under GT, since the "distance" 
between points changes under a GT when measured non-simultaneously. In this case, it is straightforward to check 
that infrared divergences do not cancel in either ED or ZA. A similar effect happens when computing the phase shift 
of the density field, which is found to be infrared divergent [67], as expected from the fact that it is not a Galilean 
invariant quantity (see Eq. (7.12)). 

We now examine the properties of the different non-linear approximation dynamics under a GT, in particular the 
momentum conservation equation, which is the one being modified from the ED in each case. For FFA, we can write 
the evolution for #(x, r) as: 

^i-?fW)=0, (7,7, 
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which is obviously non-invariant under GT. We can view the non-cancellation of infrared divergences in FFA as a 
result of the non-cancellation of the convective terms in the infrared region, since these are absent in the momentum 
equation. Therefore, the breaking of Galilean invariance in this approximation leads to infrared-divergent results. 
The case of LPA is a little more subtle. The velocity divergence equation is: 

+U(t) 0(x, r) + ^ 2 (t)Mx, r) = -V • {[v(x, r) • V]v(x, r)}, (7.18) 

the only difference with the exact dynamics being the Poisson term on the left hand side. This term is only present in 
linear theory, being replaced by zero in higher orders. Therefore, LPA effectively treats the dynamics in an "external 
field" Si(x, r) that is fixed by the initial conditions. Even though Eq. (7.18) is invariant under GT, Galilean invariance 
is "spontaneously broken" because the external force leads to a non-zero bulk velocity. By dimensional analysis, the 
induced rms flow is v rms « \\7(f>i\/H « 4>\/{%R), where R is the characteristic comoving distance between minima 
of <f>i set by the initial conditions. In fact, in LPA, particles oscillate about the minima of the linear gravitational 
potential with a characteristic velocity which can be derived for simple potentials [49]. 
In the ZA, the equation of motion for #(x, r) is: 



S0(x,r) H( 



dr 



-V-{[v(x,r)-V]v(x,r)}, (7.19) 



which is invariant under GT, consistent with the fact that loop corrections in ZA are well behaved in the infrared. 
Note that at the level of the equation of motion for v(x, r) (i.e., Eq. (2.12)), ZA is not invariant under GT; however, 
invariance is recovered at the velocity divergence level, which is what matters for this discussion. 

One can now ask whether the exact equations of motion could be modified in a different way from ZA, so as to yield 
a different non-linear approximation with infrared-convergent loop corrections. Starting from the exact dynamics, 
one can write a general class of non-linear dynamical theories of gravitational evolution in the form: 

Xi dS(^Tl + Az[v(X; r) _ v]rf(Xj r) + Ag[1 + r)]j?(Xj r) = Qj (? 2Q) 

lo "^ T) + 7i V • {[v(x, r) • V]v(x, r)} + 72 W(r) 0(x, r) + ^ 2 (r) 73 <5(x, r) = 0, (7.21) 

where 7^ (y = 0,1,2,3) and A 8 - (i = 1,2,3) are arbitrary coefficients that characterize the dynamical model. ED 
corresponds to 7^ = (1,1,1,1), ZA to 7^ = (1,1,-1/2,0) and FFA to 7^ = (1,0,-1/2,0). All of them have 
A 8 - = (1, 1, 1) (LPA cannot be cast into this form). By dimensional analysis, one could add new terms to the equations 
of motion, such as TiOS to Eq. (7.21), but their physical origin is not clear, so we do not include them. Without loss 
of generality, we can set 70 = A3 = 1. Galilean invariance then is equivalent to the conditions 

7i = l, Ai=A 2 . (7.22) 

From the generalized equations of motion (7.20), (7.21), we can derive the corresponding perturbation theory kernels 
by following the procedure of Section II B. In order for this family of theories to correctly describe linear gravitational 
instability in the limit S, 6 <C 1, we must have: 

Ai = l, 2 72 - 3 7 3 = -l. (7.23) 

Another condition comes from the fact that J(x, r) and #(x, r) vanish at the homogeneous mode, that is S(h = 
0) = #(k = 0) = 0. This plus translation invariance (which guarantees that momentum is conserved at each diagram 
vertex) implies that, e.g., ^(q, — q) = and G 2 (q, — q) = 0. Using the perturbation theory kernels for this theory, 
these two conditions lead to, respectively: 

(9 7 3 - 4A l72 )(A 2 -l) _ 

(8 72 A 1 - 1573) " ' ( ' 

7s(A 2 - 1) 
( 872Al _ i5 73 ) " °' 

which imply that A 2 = 1. This in turn guarantees the Galilean invariance of mass conservation. 
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We now calculate the divergent contribution to the one-loop variance of the density field (setting Ai = 1 but keeping 
A2 arbitrary for the moment): 



4(9 73 - 4 72 )( 7 i + 2 72 A 2 - 3 73 A 2 ) 2 
27(8 72 - 15 73 ) 2 (4 72 - 7 73 ) 



■A, 



(7.26) 



whereas for the divergent part of t^ 



we get: 



8 73 ( 7 i + 2 72 A 2 - 3 73 A 2 ) 2 
9(8 72 - 15 73 ) 2 (4 72 - 7 73 ) 



A, 



(7.27) 



(a similar calculation for the cross-correlation < 80 > to one loop does not add a new condition). Note that the poles 
of Eqs. (7.26) and (7.27), 73 = —4/3, —2 correspond to repulsive gravity and therefore are not relevant to the present 
discussion. From Eqs. (7. 22), (7. 23), (7.26), and (7.27), we conclude that under the assumptions made above, Galilean 
mvariance is a necessary and sufficient condition for the cancellation of infrared divergences in non-linear theories 
that have the correct linear limit. These equations also show the general structure of the terms that cancel each other. 
Note that all the terms in Eqs. (7.26) and (7.27) are proportional to either A 2 or 71, which are the convective terms. 
While this discussion does not constitute a rigorous proof to all orders in perturbation theory, it does demonstrate 
that Galilean invariance plays a key role in the cancellation of infrared divergences. 

We can now answer the question of whether it is possible to find a different non-linear approximation from Z A which 
is well behaved in the infrared. The restrictions imposed by Galilean invariance, translation invariance, and linear 
gravitational instability leave us with a family of non-linear theories which satisfy the constraints A 8 - = 1, 70 = 71 = 1, 
and 27 2 — 37 3 = —1. Taking 7 3 ^ means solving the Poisson equation (with a renormalized Newton's constant 
G' = 73 G). As we discussed above, this implies introducing a fundamental non-locality which greatly complicates the 
calculation of loop diagrams because of logarithmic corrections. Therefore, if 73 7^ 0, there is no reason to consider 
cases other than 7 3 = 1, which corresponds to ED. On the other hand, if we take 7 3 = 0, we arrive at the ZA. In this 
sense, ZA is the only (simple) non-linear approximation which is well-behaved in the infrared. 



We explicitly developed the general diagrammatic loop expansion for one-point cumulants of cosmological fields 
in NLCPT. We applied this formalism to calculate one-loop corrections to the variance and skewness of the density 
and velocity divergence fields in the exact dynamics of gravitational instability for scale-free initial power spectra. 
One-loop corrections for the unsmoothed density field dominate over tree-level contributions when af k, 1/2. For the 
divergence of the velocity field, this dominance does not happen until af k, 1. The results show a weak dependence on 
the power spectral index, n, induced by logarithmic terms coming from the non-locality of the perturbative solutions. 
This spectral dependence is absent at tree-level (for unsmoothed fields) , because tree diagrams correspond to averaging 
out the tidal field which contains the non-local contribution. In fact, the dynamics at the tree-level has been shown 
to be equivalent to the spherical collapse model, and can be characterized by the "monopole" (angular average) of 
the perturbation theory kernels [7,17,18]. 

We also calculated loop corrections in different non-linear approximation schemes: the Zel'dovich approximation 
(ZA), linear potential approximation (LPA) and the frozen-flow approximation (FFA). Of these three theories, ZA 
shows the closest loop corrections to the ED, followed by LPA, and then FFA, in qualitative agreement with the 
known results at the tree-level. Calculations in this set of approximations are much simpler, due to the absence 
of non-local contributions which necessarily appear when solving the Poisson equation beyond linear perturbation 
theory. This advantage of course dissapears when one considers the calculation of higher-point cumulants, since now 
the statistic itself becomes non-local. Of all the non-linear approximations, ZA is the only one which gives loop 
corrections convergent in the infrared; we showed that this fact is related to the invariance of the equations of motion 
under a Galilean transformation. 

These results show that loop corrections are generally not negligible in the weakly non-linear regime a < 1 for 
unsmoothed fields. However, to see whether NLCPT can describe the a dependence of one-point cumulants, the 
effects of smoothing must be included in order to compare the predictions with numerical simulations. Results 
obtained recently for one-loop corrections to the variance with gaussian smoothing show very good agreement with 
N-Body simulations [32]. This is certainly encouraging news, and gives support to the perturbative approach to 
gravitational instability. In a subsequent paper, we will present loop corrections for smoothed cumulants and multi- 
point correlation functions. 



VIII. SUMMARY AND CONCLUSIONS 
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APPENDIX A: ANGULAR INTEGRALS 



In this appendix we calculate the necessary angular integrals for loop corrections to 1-point cumulants. For correc- 
tions of order O(af) to cumulants in the non-linear approximations we need the integral: 

J dfii J di2 2 J di2 3 ( qi • q 2 ) a (q 2 • q 3 ) 6 (q 3 • qi) C , (Al) 

where a, b and c are non-negative integers. This integral, and its generalizations to higher number of wave-vectors 
can be calculated by using the result that: 

,„ I q 1 ' ' Ai. i if v is even, , . 

q Hn 1 ' p 10 if p is odd, v ; 



where: 
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A, 



iy J^^t(»'i)t(»2) • • A(ip-i)T(i P )> (A3) 



and Sij denotes the Kronecker symbol. These results follow from antisymmetry under inversion of a component of q, 
symmetry under permutations of the i's, and rotational invariance. The coefficient in Eq. (A3) is obtained from the 
evaluation of a particular component of Eq. (A2). Eqs. (A2) and (A3) are the only integrations needed in the case of 
the non-linear approximations considered in the main text. Note that the result of the angular integrations gives a 
function in which the q's are decoupled, so the radial integrations reduce to a product of one dimensional integrals. 

Loop corrections in ED require the calculation of more complex angular integrals, due to the presence of angular 
variables in denominators coming from the Poisson Green's function. This introduces the complication that the 
integration does not decouple into one-dimensional integrals, a reflection of the non-locality introduced by the exact 
dynamics of gravity. (This results in the spectral index dependence of loop corrections.) For corrections up to order 
0(af) we need: 

h ^[dnJ dn 2 S d n 3 (qi-q 2 ) a (^«i3) t (q3-qi) e (A4) 

J J J (qi + q2r 

We choose the polar axis in the direction of qi and the «-axis so that q 2 lies in the x-z plane. Then: 

qi = <Zi(0,0,l) , (A5a) 
q 2 = ?2(sin 9, 0, cos 9) , (A5b) 
q 3 = 53 (sin u> cos <f>, sinw sin <f>, cos u>) , (A5c) 

and therefore: 

fi fi [-2tt x a \\J \ — x 2 \J \ — y 1 cos <f> + xy I y c 

h = ^ 2 ql +a q a + b q b + c dx dy dj> —± ^— 3— >- , (A6) 

J-l J-i Jo q( + q$ + 2q 1 q 2 x 

where x = cos 9 and y = cos w. Since b is non-negative, we can use the binomial expansion to get: 

/i = 87rV«2 a+6 «S +e E(-) / dx „* + „* + Z„r / d yy c+b ~ l ( l -y 2 y /2 d^^y. (at) 

i=0 Vv J -1 9i ~r ( h t L( l\'n x J -1 Jo 

Now note that for i = 2k: 
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(cos.) 2 * = 2^"^ 



(2A)!! ' 

the corresponding integral being zero when i is odd. This result allows further binomial expansion in (A7) to get 



(A8) 



a+b-2i 



dx 



qf + qi + 2qiq 2 x 



dy y 



,c+b-2j 



(A9) 



where [6/2] denotes the integer part of 6/2. The integral over y is elementary, and using Eq. (B3) (see Appendix B) 
we arrive at the desired final expression: 



"=^i(i)¥s:(!)^') 

k = V 7 y ' i = V 7 j = yJ/ 
q'{ + qly+b-'ii-'ir-i | ( q l + q ?y+b-2i \ qi + q 2 \ 



-iy+j 



(b + c-2j + 1) 



[(a+6-2i-l)/2] 



E 



i 



2r + 1 



2gi^2 ^ ' ^ 2gig 2 ^ \qi - q-i | 

if (6 + c) is even, and zero otherwise. Another angular integral of interest is: 

(qi • q2) a (q2 • qs) 6 (q3 • qi) 



h = / rffii / drt 2 / dn 



(qi + q2) 2 (qi + qs) 2 



(A10) 



(All) 



Following the same steps that lead to Eq. (A10) we obtain: 

ft) £(;)<-: 



c+a-2 a+b-1 b+c-1 

1i 'h </3 



1 



2r+l V 2 qi q 2 



q\ + ql\ a+b ~ 2i Xn l?i + ?2| 



ki - 92 1 



x * 7 — V ? 

[(6+c-2j-l)/2] 



i+j 



[(a+6-2i-l)/2] 

E 

r = 



E 

s = 



9 2 + g|\ 6+c-2j-2 S -l 1 

2 qi q 3 7 2s + 1 



l?i - 93 1 

We also need angular integrals involving coupling of three wave- vectors in the denominator: 

(qi -q2) a (q2 -q3) 6 (q3 • qi) c 



h = / rf^i / ^2 / dn 



(qi + q2 + q3j 



q\ + g2^a+6-2i-2r-l 
2^192 



9i + Qs\ b+C ~ 2j 
2gig3 



(A12) 



(A13) 



In this case it is more convenient to choose the x-z plane as the plane spanned by qi and q 2 , and take the z-axis to 
be along Q = qi + q 2 . With these conventions we have: 



qi = gi(sin 6*i , 0, cos 6\) , 

q 2 = g 2 (-sin6> 2 ,0,cos6> 2 ) , 

q3 = g3 (sin 6s cos <f>, sin 6s sin <f>, cos 6s) 



(A14a) 
(A14b) 
(A14c) 



All the necessary dot products can be expressed in terms of <f>, z = cos 6s and x = cos 6 where 6 is the angle between 
qi and q 2 : 



qi • q 2 = qiq 2 x 



q 2 • q 3 = q 2 qs[~ qi V 1 - x 2 V 1 - z 2 cos <j) + (q 2 + q 1 x)z]/Q 
q 3 • qi = q3qi[q2\/l ~ x 2 \/l - z 2 cos <j> + (q 1 + q 2 x)z]/Q, 



(A15a) 
(A15b) 
(A15c) 
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where Q 2 = q 2 + q 2 + 2q\q 2 x. We therefore have: 

-,; 

qw/l — x 2 \J\ — z 2 cos <f> + (q 2 + 



.•1 /-l /-27T X" 

2c+aa+66+c 



x [?2\/l - * 2 Vl - ^ 2 cosc?i + (gi + g 2 «)^] • (A16) 
Using binomial expansions, Eq. (A8) and the following identity for any function F: 

be b [(«'+c)/2] 

E E ^') = E E F{i,2k-i), (A17) 

i = i = i = fe = [(i + l)/2] 

we get after some lengthy algebra: 

i = fe = [(i + l)/2] v 7 v 7 v ; t = u=0u=0 \ / \ / \ 



v c+ti-2fe-t + 2i 6-ti+t + 2fe-2i 
x< /l 92 



1 j.a+t+u + 2); * /£\ /•! 



r = 

Now, using Eq. (B3) we have: 

b [(i+c)/2] c-2fe+i 6-i k /JA / \ /; \ /l -\ / o; i -\ /; \ /o; mi 

/s = i6A^2 + ^ +c E E E EEE c ^ m A " A (c ~ Am ( } 



x (_ 1 ^+«+r g c +tl -^- t + ^^o 



jj\2k-ij\vj\ u \ t J\rJ (2k)\\ 

! = k = [(i + l)/2] t = u = 0v = 0r = w v /\/\ /\ / \ / \ J 

r [(b+c-2k + 2r-l)/2] , .b+c+2r-2k-2s-l ^ b+c-2k + 2r-2s-l 

^ \2q^J 2s + 1 ^ 

s = v " / w = 



i+v+r c+u-2k-t + 2i b-u+t + 2k-2i 



h + c -2k + 2r-2s-l\ 2w _ x f 1 J x °+t+«+^ s^b+c+^^k b+c-2k+2r , h + c _ 2k + 2 . ^ 



w = 



1 dx x «+t+«+2« \Q + q 3 \ 
Finally, noting that: 



f 1 dxx™ \Q + q 3 \ 2 A/mV lvn -A, 2, 2^ m - A / |91+921 , 2(A _ C) |z + g 3 | 



(A19) 



(A20) 



we arrive at the final expression 

[(i+c)/2] c-2fe+i b-i k k / ,n 

JJ\2k-iJ\vJ\ u )\ t J\rJ (2k)\\ 



8 = fe = [(i + l)/2] t = u=0»=0r=0 

[(6+c-2fe + 2r-l)/2] / 1 N b+c+2r - 2k - 2s - 1 1 6+c- 2fe + 2r- 2s - 1 

w / i +r „c+ti-2fe-t + 2j' 6-ti+t + 2fe-2i 



/ ^ \ u-\-c-\- zr — zk — zs — l ^ 

^ \2q^J 2s + 1 ^ 

/6 + c-2Ar + 2r-2s-l\ q 2 ™' 1 , q 2 + ql 

' ' ^ 3 M(a + t + n + 2t;, k - r + s + w + 1, ^ 2 ;-!,!) 



V w / {2q 1 q 2 ) k - r + s + w + 1 y ' ' 2 gi g 2 

+ U 3 J ^ V » 7 (2q 1 q2) a + t + u + 2 ^ 1 ^ I /» 

x( _ 1)a + t+U + 2„-A (g 2 + q 2 2) a +t+ u + 2v-h AJ{2h _ 2k + 2r _ 2w j ^ ; | gi _ g2 |, | gi + g2 |) 



(A21) 
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where we have used the definitions: 



M(m, c, a ; a,b) = / dx 



(x + a) c ' 

with M(m, I, a ;a,b) = M(m, I, a ; b) — M(m, I, a ; a) (see Eq. (B3)), and: 



(A22) 



J(m,a;a,b)= / cfe « m In \x + a|, 



(A23) 



with AJ(m, a ; a, 6) = J(m, a ;a,b) — J(m, —a ; a, 6). These integrals will be calculated in the next appendix since 
they are also needed for radial integrations. 

From Eq. (A21) one can get the result for the angular integral: 



h = / rffii / dtt 2 / dtt 



(qi + q2 + q3) 2 (qi + q2) 2 

since the only difference with Eq. (A13) is the additional factor of l/Q 2 . We get: 

6 [(i+c)/2] c-2k+i b-i k k 

Si \2k — i ) \v 



(A24) 



,4 *Vwe "f T' £Ei:C)( 2 ;_,)C)Cr)C" 2 < +, )C) W 

8 = fe = [(i + l)/2] t = u=0u=0r=0 v/v /\/\ / \ / \ / \ I 



-\\i+v+r c+u-2k-t + 2i b-u+t + 2k-2i 



[(6+c-2fe + 2r-l)/2] . ^ . 6+c+2r- 2fe - 2s - 1 ^ 6+c- 2fc + 2r- 2s - 1 

^ V2glJ 2s + 1 E 

s = v " / iu = 



b + c-2k + 2r -2s - 1 

w / (2gig 2 



A; — r + s+u/ + 2 



M(a + t + n + 2v, k - r + s + w + 2, Ql . + ?2 ;-!,!) 



6+c+2r _ 2fc 6+C -2* + 2r ^ + £ _ ^ + ^ 



E 



2 gf" 1 



(2gig2 



ia+t+u + 2u + 



a+i+u + 2i/ 

r E 



2gi^2 

a + t + m + 2i> 
h 



W = V " ' V a-^a^/ />_(] 

;( _ 1)a +t+ U + 2„-A (g 2 + q 2 )a+ t + u + 2v-h AJ{2h _ 2k + 2r _ 2w _ 2jq3 . \ qi _ g2 | ; \ qi + q2 \ 



(A25) 



APPENDIX B: RADIAL INTEGRALS 



We now proceed to calculate the radial integrals involved in loop corrections to cumulants in the ED case. As the 
results from the angular integrations show, most of the radial integrations needed involve logarithms. The number of 
radial integrals involved is quite large, especially because in general radial integrals do not decouple into products of 
independent one-dimensional integrations. However, all these integrals can be calculated in terms of the ones given 
in this appendix using integration by parts and partial fraction decompositions. Many of the integrals of interest are 
generated by: 

f x m (-l) c_1 9 C_1 

M(m, c,a ;x)= dx - — ■ — = — M(m, I, a ;x). (Bl) 

J (x + a) c (c— 1)1 oa c ~ L 

The integral for c = 1 can be solved by noting the following identities for m > and a ^ 0: 

m 

Y,x i - 1 {-a) m - i , (B2a) 



x"" {-a 



x + a x + a 

i=i 



1 1 1 

y ■, ^ • • (B2b) 



x' n (x + a) (— a)' n (x + a) ^—^x'[—a 



Therefore 
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_^ (—01 Y J X J 

M(m,l,a ;x) = 2_^± '— h (-a) m In \x + a\ , m>0, (B3a) 



/ _ a) m+j 



M{m,l,a ;x) = ^ 1 , ' . — + {-a) m ( In \x + a\ - In \x\j ,m<-l. (B3b) 

By taking derivatives with respect to a (denoted by d a ) according to Eq. (Bl) one can generate M(m, c,a ;x) with 
c ^ 1 . We get : 

M(m, c,a;x) = ^2 * 1J 0*" V" 1 ) + ( / n , 1 («) m In |z + «l + £ • T^YT 

^ (m - j + l)(c- 1)! (c-lj! J /(c-!)! 

x l J 9°- 1 - J "(a m ) ,m>0, (B4a) 

M(m,c,a;x)= ]P \ j 5°~ 1 (q~ j ) + , n , 1 («) m ( In \x + «| - In \x\) + ]P . 

^ (m + j)(c-l)! (c-lj! V ) j J 

* (c-l)! x + a ' m ~ ( B4b ) 

Another important integral is: 

J(m, a ; x) = dx x m In \x + a\, (B5) 



which can be solved by shifting the integration variable x — > x + a and then using integration by parts for m > 
and m = —2. The m < —2 integrals are generated by derivatives of J(m, a ; x — a) with respect to a. The m = —1 
case deserves special treatment and needs the introduction of a new function, the dilogarithm. Shifting variables back 
we get: 



(B6a) 



m + l _ l_ a \m + l ™ / \ / , a \j + l 

J(m, a;x) = [ -p In \x + a\ - V (-a) m "-> V "\ 9 , m > , 

1 ; m + i 1 1 ^Vj/ (J + 1) 2 " 

</(— 1, a ; = In \a\ In |a;| — Re[Li 2 (— x/a)\ , (B6b) 

x m + l _ (_ a \m + l (_ a )m + l (_ a )m + l "™^ 2 (_ a )j 

JSm, a ; x) = In \x + a 1 + ^ '- InM-- '- > ^ J t- ,m<-2, (B6cj 

v; m+l 1 1 m+l 11 m+l ^ > _ > v ; 

where we have used Re[ln(a;)] = In \x\, and Li2 denotes the dilogarithm, defined by [69]: 

U 2 (x) = - T dz ln(1 ~ z) . (B7) 



■jo <■ 

This definition is extended to higher order polylogarithms in the following way [69]: 

Li n (z) = f dz Lln " l(z) , (B8) 
Jo z 

where Lii(a;) = — ln(l — z)/z. Polylogarithms have the series expansion for small argument (x < 1): 

oo ^ 

Ll nW=E^' ( B9 ) 



which gives the connection to the Riemann zeta function, Li n (l) = C( n )- Polylogarithms of real argument x are 
real as long as x < 1, and for x > 1 we have Im[Li n (a;)] = — irlxi" -1 (x)/ (n — 1)!. To carry out series expansions of 
polylogarithms in powers of e/k c we also need the so-called functional equations [69]: 
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Li 2 (a; 
Li 2 (-x 

Li 3 (x 
Li 3 (-x 

\A^(x 
\A/±( — X 



-Li 2 (l/x) 



■ In 2 Sx) 



7T 

y 



in In ; 



, x > 1, 



-Li 2 {-l/x) - ^\n{x) - y, 



Li 3 (l/x) 



ln 3 (x) 



ln(x) 



1 1 2 
— 7T 111 X 



, X > 1, 



1 3 7T 2 

Li 3 (-l/a;) - - In (a;) - y In (a;) 



-Li^l/x) 



1 

24 



ln 4 (a;) 



■ In 2 (a;) 



6 



* « , 3 



, a; > 1, 



1 7T 2 77T 4 

-Li 4 (-l/a;) - — ln 4 (a;) + y ln 2 (a;) - — , 



(BIO) 
(Bll) 
(B12) 
(B13) 
(B14) 
(B15) 



which show the behavior of polylogarithms for large positive and negative arguments. The integrals involving poly- 
logarithms combined with powers and/or logarithms that appear in the calculation of 1-loop corrections to a 2 and S3 
can be obtained in terms of M , J and the following last "master" integral given by: 

T(m,a,/3 ,x) = J dx x m \n\x + a\ \n\x + ,B1«> 

which can be easily solved in terms of J integrals by integrating by parts when In this last case, unless a 

or [3 are zero, the integration by parts does not work and one has to calculate it explicitly. By factoring out a and [3 
from inside the logarithms one finds: 



T(— 1, a, [3 ; x) = In \a\ In | 



In \x\ - Re[Li 2 (-^)] In |/?| - Re[Li 2 (-^)] In 



~ xla dy 

— In 1 1 — y I In |1 - cy\ 



y 



(B17) 



where c = a//3. This last integral is solved by using that 

l — cx 



— ln(l - x) ln(l 

x 



Li 3 



1 — x 



Li 3 (l/c) + Li3(l)-Li3(l-a0-Li3 



2 

+ ln(l - cx) [Li 2 (l/c) - Li 2 (a;)] + ln(l - x) [Li 2 (l - cx) - Li 2 (l/c) + y] . (B18) 



6 

Finally, we must mention that special cases arise when the constants a or [3 become zero or such that there is a 
singularity in the region of integration. Most of the time these situations can be solved by taking the appropriate 
limits (by using the functional equations for the polylogarithms, for example). There are cases, however, where 
individual integrals are singular, but these divergences can be regulated by introducing another small parameter, say 
S, which gets canceled at the end of the calculation by taking S e <C k c . In practice, the symbolic integration 
routine developed checks for singularities in the integration region and introduces infinitesimal cutoffs accordingly. 
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